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Abstract 

We review the current understanding of the diffuse gamma-ray background (DGRB). 
The DGRB is what remains of the total measured gamma-ray emission after the sub¬ 
traction of the resolved sources and of the diffuse Galactic foregrounds. R is interpreted 
as the cumulative emission of sources that are not bright enough to be detected individ¬ 
ually. Yet, its exact composition remains unveiled. Well-established astrophysical source 
populations (e.g. blazars, misaligned AGNs, star-forming galaxies and millisecond pul¬ 
sars) all represent guaranteed contributors to the DGRB. More exotic scenarios, such as 
dark matter annihilation or decay, may contribute as well. In this review, we describe 
how these components have been modeled in the literature and how the DGRB can be 
used to provide valuable information on each of them. We summarize the observational 
information currently available on the DGRB, paying particnlar attention to the most 
recent measurement of its intensity energy spectrum by the Fermi LAT Collaboration. 
We also discuss the novel analyses of the auto-correlation angular power spectrum of 
the DGRB and of its cross-correlation with tracers of the large-scale structure of the 
Universe. New data sets already (or soon) available are expected to provide further 
insight on the nature of this emission. By summarizing where we stand on the current 
knowledge of the DGRB, this review is intended both as a useful reference for those 
interested in the topic and as a means to trigger new ideas for further research. 
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1. Introduction 

The first full-sky image of gamma-ray emission was obtained in 1972 by the OSO-3 
satellite. It consisted of 621 events detected above 50 MeV [1]. Since then, telescopes 
with lower sensitivities and better angular and energy resolutions have signihcantly 
improved our understanding of the gamma-ray Universe. The all-sky maps produced 
by the Fermi Large Area Telescope^ {Fermi LAT from now on) after more than 6 years 
of data taking contain more than 5 million events above 1 GeV. These maps exhibit a 
rich morphology: along the Galactic plane, the diffuse Galactic foreground is the most 
evident feature and, overall, it accounts for ~ 80% of the detected gamma rays. This 
diffuse radiation is produced by the interaction of cosmic rays (CRs) with the Galactic 
interstellar radiation field and with the nuclei of the Galactic interstellar medium. Also, 


^ http://fermi.gsfc.nasa.gov/ 
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the Third Fermi LAT catalog (3FGL) reported the detection of 3033 sources throughout 
the sky [2], Extended gamma-ray emitters are discussed in Ref. [3]. Other structures, 
e.g. the Fermi bubbles [4], represent more complex phenomena, whose emission is not 
fully understood yet. 

In order to reproduce the data from the Fermi LAT it is necessary to include one 
additional contribution, i.e. a diffuse and nearly isotropic emission called the Diffuse 
Gamma-Ray Background (DGRB) [1, 5, 6, 7, 8, 9]. The DGRB is thought to be pre¬ 
dominantly of extragalactic origin; gamma-ray sources with a flux smaller than the 
sensitivity of Fermi LAT are not detected individually, producing instead a cumula¬ 
tive diffuse glow that contributes to the DGRB. Unresolved blazars [10, 11, 12, 13, 14, 
15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], misaligned Active Galactic Nuclei (MAGNs) 
[26, 27, 28, 29], star-forming galaxies (SFGs) [30, 31, 32, 33, 34, 35] and millisecond 
pulsars (MSPs) [36, 37, 38] are guaranteed components to the DGRB. More uncertain 
source classes, e.g. galaxy clusters [39] or Type la supernovae [40, 41], may also play a 
role, together blue with diffuse phenomena, e.g. the radiation produced by the interac¬ 
tion of Ultra-High-Energy Gosmic Rays (UHEGRs) with the Extragalactic Background 
Light (EBL) [42, 43].^ It is possible to estimate how much different source classes con¬ 
tribute to the DGRB, but its exact composition remains one of the main unanswered 
questions of gamma-ray astrophysics. Finding a definitive answer would constrain the 
faint end of the luminosity function of the DGRB contributors. Indeed, the study of the 
DGRB may represent the only source of information about those objects that are too 
faint to be detected individually. 

The DGRB may also shed some light on exotic Physics as, e.g., on the nature of Dark 
Matter (DM): a huge experimental effort is currently devoted to the so-called indirect 
detection of DM, i.e. the search for particles (e.g. gamma rays, neutrinos, positrons 
or anti-protons) produced by the annihilations or decays of DM. Such a signal would 
constitute the first evidence that DM can interact non-gravitationally and it would 
represent an enormous step forward in our understanding of its nature [44, 45, 46, 47]. 
Targets like the center of the Milky Way (MW), local satellite galaxies or nearby galaxy 
clusters are considered optimal, thanks to the intensity of the expected DM signal and/or 
to the absence of significant competing backgrounds. Yet, no signal has been robustly 
associated with DM. Then, if DM annihilates or decays producing gamma rays and 
such a signal has not been detected up to now, it is most probably unresolved and it 
contributes to the DGRB. Looking for the features of a DM component in the DGRB, 
one can, then, hope to finally unravel the long-standing mystery of the nature of DM 
[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 45, 60, 61, 62, 63, 64, 65]. 

The most recent measurement of the intensity of the DGRB has been performed 
by the Fermi LAT, in the range between 100 MeV and 820 GeV [9] (see also Refs. [1, 
5, 6, 7, 8] for the previous measurements). Valuable information on the nature of the 
DGRB can be extracted from its intensity energy spectrum. Its steepness could indicate 
whether a particular class of sources dominates the emission. Moreover, the transition 
between two energy regimes dominated by different classes could, in principle, give rise 


^Different names have been used in the literature to denote the DGRB, e.g. Extragalactic Gamma-Ray 
Background or Isotropic Gamma-Ray Background. We believe the denomination used in this review is 
more precise since, as we will see in the following sections, the DGRB may be not entirely extragalactic 
and since it exhibits a certain amount of anisotropy. We also note that, in Ref. [9], the Fermi LAT 
Collaboration uses the name Extragalactic Gamma-Ray Background to characterize the cumulative 
emission of all sources (both resolved and unresolved), while Isotropic Gamma-Ray Background refers 
to the unresolved component only, i.e. what we call DGRB (see Sec. 2). 
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to breaks and features in the energy spectrum. Yet, since the intensity of the DGRB is 
only sensitive to the sum of its contributions, ultimately there is only a limited amount 
of information that can be extracted from it. 

Fortunately, due to its excellent sensitivity and exemplary angular resolution, the 
Fermi LAT marked the beginning of an era in which the intensity energy spectrum is 
no longer the only observational data available for the study of the DGRB. In 2012, 
the hrst measurement of anisotropies in the DGRB was reported [66], and the detection 
of a non-null auto-correlation angular power spectrum (APS) provided complementary 
constraints on the composition of the DGRB [67, 68, 63]. Moreover, in Ref. [69], the 
authors used the photon-count probability distribution measured after 11 months of 
Fermi LAT data to constrain the contribution of unresolved blazars to the DGRB. More 
recently, the cross-correlation of the DGRB anisotropies with observables tracing the 
Large Scale Structure (LSS) of the Universe has also been considered. In Ref. [70] the 
authors measured the 2-point correlation of the DGRB with 5 different galaxy catalogs 
and they reported a signal with 4 of them. Ref. [71], instead, cross-correlated the DGRB 
with the cosmic shear observed by the Ganada France Hawaii Telescope Tensing Survey 
(CFHTLenS) and found no significant detection. These works, together with Refs. [72, 
73, 74, 75] , have proved that the study of the cross-correlation of the DGRB with the LSS 
is a very powerful strategy that may provide access to components of the DGRB that are 
only subdominant in the intensity energy spectrum or in the auto-correlation APS. In 
particular, the technique has the potential to deliver the hrst detection of DM-induced 
gamma-ray emission. 

In the near future, the data on the DGRB is expected to further increase, due 
to the extended run of the Fermi LAT and to the fore-coming Gherenkov Telescope 
Array (GTA) [76, 77]. Other frequencies and other messengers will also be crucial to 
improve our modeling of the DGRB and to extract complementary information about 
its nature. The scenario is rapidly evolving and the study of the DGRB is quickly 
becoming a standard tool for the characterization of unresolved astrophysical sources 
and of a potential DM-induced gamma-ray signal. 

Therefore, we believe that this is the right moment to summarize where we stand 
in our understanding of the DGRB. In this article we will survey the data available 
on the DGRB at present and we will discuss how these observations have been used 
to constrain the nature of the emission. We will also enumerate the classes of sources 
or emission mechanisms that have been proposed as contributors to the DGRB. By 
sketching a snapshot of the state-of-art on the DGRB circa 2015, we intend to provide 
the community with a reference point from which to build on. 

We end by noting that the DGRB is intrinsically an analysis- and time-dependent 
quantity. Indeed, its intensity depends on the sensitivity of the telescope employed to 
detect it and on its instrumental capability to resolve sources. Even with the same detec¬ 
tor, an increase in statistics or, in general, any improvement in the detection sensitivity 
will result in a different DGRB. In the following sections, every time we mention the 
DGRB we will make sure to specify which measurement of the DGRB we refer to. 

The paper is organized as follows: in Sec. 2 we focus on the intensity of the DGRB. 
Sec. 2.1 reviews the recent measurement of the DGRB energy spectrum by the Fermi 
LAT, while Secs. 2.2 and 2.3 are devoted to the description of the sources that have 
been proposed as contributors to the emission: astrophysical objects are studied in Sec. 
2.2, while Sec. 2.3 discusses the case of the DM-induced emission. Sec. 3 reviews the 
Fermi LAT measurement of the auto-correlation APS of anisotropies and its impact 
on our understanding of the DGRB. The topic of Sec. 4 is the measurement and the 
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interpretation of the photon count probability distribution, while in Sec. 5 we investigate 
the cross-correlation with probes of LSS, namely galaxy catalogs in Sec. 5.1, cosmic 
shear in Sec. 5.2 and other observables in Sec. 5.3. Finally, we present our conclusions 
in Sec. 6. 

2. The intensity energy spectrum 

The subject of this section is the intensity energy spectrum of the DGRB. The first 
subsection (Sec. 2.1) summarizes the most recent measurement of this observable, i.e. 
the one performed by the Fermi LAT in Ref. [9]. The following subsections (Sec. 2.2 and 
Sec. 2.3) present the different source classes and emission mechanisms that have been 
proposed to interpret the observed emission. For each population we summarize how the 
DGRB energy spectrum has been used to learn about the properties of the gamma-ray 
emitters. 

2.1. The new Fermi LAT measurement of the Diffuse Gamma-Ray Background 

The Fermi LAT [78] on board the NASA Fermi satellite started its scientific oper¬ 
ation on August 2008 and, since then, it has revolutionized our knowledge of the most 
violent phenomena in the Universe. The Fermi LAT covers 4 decades in energy, from 
few dozens of MeV up to the TeV regime. It has an unprecedented sensitivity (~ 30 
times better than its predecessor EGRET) and an extremely large field of view reaching 
almost one fifth of the sky. Its angular resolution is of about 0.8° at 1 GeV and better 
than 0.2° above 10 GeV.^ Such superb capabilities allowed the discovery of hundreds 
of new gamma-ray sources and an impressive cartography of the Galactic GR-induced 
gamma-ray diffuse emission that reaches, for the first time, energies greater than 10 GeV 
[79]. 

In addition to these achievements, the analysis of 10 months of data delivered the 
first Fermi LAT measurement of the DGRB energy spectrum at Galactic latitudes, b, 
greater than 10 degrees [8]. Such a measurement was performed between 200 MeV and 
102 GeV and it represents the third independent observation of the DGRB, after the one 
by the SAS-2 satellite in 1978 between 40 and 300 MeV [80] and that by EGRET between 
40 MeV and 10 GeV in 1998 [6]. The energy spectrum of the DGRB as measured by the 
Fermi LAT in Ref. [8] is featureless and it can be well described by a single power law 
with a spectral index of 2.41 ± 0.05. This is significantly softer than both the DGRB 
initially reported by EGRET [6] and the revised estimate of Ref. [81] from the same 
data set."^ 

More recently, a new measurement of the DGRB energy spectrum has been performed 
by the Fermi LAT [9]. The analysis used 50 months of data with \b\ > 20°, it employed a 
dedicated event selection and it took advantage of improvements in the determination of 
the GR background and of the diffuse Galactic foreground. The measurement, denoted 
by red data points in Fig. 1, now extends down to 100 MeV and up to 820 GeV. It 
represents the most complete and accurate picture that we currently possess of the 
DGRB intensity. Interestingly, the DGRB now exhibits a high-energy exponential cut¬ 
off at 279±52 GeV (for the baseline model of Galactic diffuse foreground used in Ref. [9]), 
and it is well described by a single power law with a spectral index of 2.31 ±0.02 at lower 


®These values refer to the 68% containment radius. 

^However, note that, in the re-analysis of Ref. [81], the energy spectrum is well described by a single 
power law only up to 2 GeV. 
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energies. The cut-off is compatible with the attenuation expected from the interaction of 
high-energy photons with the EBL [82, 83, 84, 85, 86, 87, 88] over cosmological distances 
[25]. The largest systematic uncertainty (represented by the red shaded region in Fig. 1) 
ranges between a factor of ~ 15% and 30% (depending on the energy range considered) 
and it comes from the modeling of the Galactic diffuse emission. 

For the purposes of this review, it is convenient to briefly summarize here the main 
steps followed in Ref. [9] to measure the DGRB. Two different sets of selection cuts are 
applied to the gamma-ray data, optimized over different energy ranges, namely below 
and above 12.8 GeV. Different selection criteria are required because of the energy de¬ 
pendence in the composition of the CR-induced backgrounds and the way GRs interact 
with the detector. For the low-energy data sample, the all-sky Galactic diffuse emission 
is modeled as a sum of templates obtained with GALPRDP^ and the point-like sources 
are modeled following the information in the Second Fermi LAT catalog (2FGL). Addi¬ 
tional templates are used to model subdominant contributions from the Loop I large-scale 
Galactic structure and from electrons interacting with the Solar radiation field through 
Inverse Compton (IC). A fully isotropic template is also included, in order to describe 
the DGRB and the residual CR contamination. The fit to the data determines the nor¬ 
malizations of the different templates and it is performed independently in each energy 
bin considered. At the highest energies, where the statistic is scarce, the normaliza¬ 
tions of the templates for the Galactic diffuse emission are fixed to the best-fit values 
obtained at intermediate energies (between 6.4 and 51.2 GeV). The spectral shapes of 
their emission are also fixed to those predicted by GALPROP. Above 12.8 GeV, then, the 
normalizations of the point-like sources and of the isotropic template are the only free 
parameters in the fit. 

Monte Carlo simulations are used to estimate the residual CR contamination. The 
energy spectrum of the DGRB is, finally, obtained by subtracting the CR contamination 
from the isotropic component determined in the template fitting. The systematic error 
induced by the imperfect knowledge of the Galactic diffuse emission is estimated by 
repeating the whole procedure for three benchmark models of Galactic foreground and 
for different values of the parameters controlling the propagation of GRs in the MW (see 
Ref. [9] for further details). 

2.2. The astrophysical components of the Diffuse Gamma-Ray Background 

In this section we describe the classes of astrophysical sources and emission mech¬ 
anisms that have been proposed as contributors to the DGRB over the years. Well- 
established astrophysical populations, whose brightest members have been robustly 
detected, represent “guaranteed” components to the DGRB. They are discussed in 
Sec. 2.2.1, Sec. 2.2.2, Sec. 2.2.3 and Sec. 2.2.4, which are devoted, respectively, to blazars, 
MAGNs, SFGs and MSPs. Then, in Sec. 2.2.5, we turn to more speculative scenarios. 
We do not discuss the possibility of a DM-induced contribution to the DGRB, since that 
is the subject of the following section (Sec. 2.3). 

We start by presenting a formalism that will be adopted throughout the manuscript 
for the description of a generic population of sources. The sources are supposed to be 
characterized by a measurable quantity Y. In most cases Y will be the blue gamma-ray 
luminosity of the source but, in some instances, it will indicate another parameter 
such as, e.g., the mass of the galaxy or of the DM halo hosting the gamma-ray emitter. 
The differential gamma-ray flux d^/dEdTl (i.e. the number of photons per unit area. 


CR propagation code available from littp://galprop.stanford.edu. See also Ref. [79]. 
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Figure 1: In red, the DGRB intensity energy spectrum, labeled “Isotropic Gamma-Ray Background 
(IGRB)” in the figure. Model A from Ref. [9] is used for the diffuse Galactic foreground. The spectrum 
is compared to the previous Fermi LAT measurement of the DGRB from Ref. [8] (blue bands) and to 
the total “Extragalactic Gamma-ray Background (EGB)” (black data points). The latter is defined here 
as the sum of the DGRB and of the resolved sources at |6| > 20° (shown in gray). The comparison 
shows that, above 100 GeV, ~ 50% of the Extragalactic Gamma-ray Background is now resolved into 
individual LAT sources. Yellow and red shaded bands represent the systematic error associated with 
the uncertainties in the modeling of the Galactic diffuse emission. Taken from Ref. [9]. 


time, energy and solid angle) expected from the unresolved objects in such a population 
can be written as follows: 


dEdQ 


(Eo) 


r^min r^max C 

I ■""'ir 

-Zmax ^ min h 

/ _ np,r,r) \ 

\ ^max / 


dr „ (z, y,r) FEoiz, y, r) x 


dVdYdr 
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where 2 ; is the redshift and T is a parameter that characterizes the shape of the energy 
spectrum of the sources. The domain of integration will depend on the particular pop¬ 
ulation considered. The factor dN/dVdYdV indicates the comoving number density of 
sources per unit Y and T, while dV/dzdVt is the comoving volume per unit redshift and 
solid angle. The factor Feq{z^Y^T) is the gamma-ray flux (at energy Fq) produced by 
the source identified by the value T, located at redshift 2 ; and with an energy spectrum 
characterized by T. If sources are described by their gamma-ray luminosity {Y = L^), 
Feo can be written as follows: 


Eeo 


(1 -|- z)‘^Ly 

AirDLizY 


dEoCr), 


( 2 ) 


where Dl(z) is the luminosity distance at redshift 2 ; and gEo(^) parametrizes the gamma- 
ray energy spectrum. For other choices of Y, the relation between Y and Feq needs to 
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be specified case per case. 

Thus, the first line of Eq. (1) indicates the cumulative emission expected from all 
the sources located between Zmm and Zmax and with characteristics between (Tmin, T min ) 
and ( Tma x. Tmax)- The quantity T, T) is called sky coverage and it is defined as the 
area O ma x surveyed by the telescope multiplied by the detection efficiency. It represents 
the probability for a source characterized by parameters {z,Y,T) to be detected and it 
encodes the sensitivity of the instrument, i.e. 11 ( 2 ;, y,r) is equal to zero if the source is 
too faint to be detected. It also accounts for any potential selection effect. An estimate of 
i}{z, Y, r) for the high-latitude blazars in the first Fermi LAT catalog (IFGL) is provided 
in Ref. [18]. Then, the factor (I — 11(2;, T, r)/nmax) in Eq. (1) has the effect of selecting 
only the sources that remain undetected. Finally, the exponential accounts 

for the effect of the EBL attenuation. 

It is also convenient to define a generic expression for the differential source count 
distribution dN/dS, i.e. the number of sources per unit solid angle and per unit flux S: 



( 3 ) 


where S is the gamma-ray flux associated with the source characterized by {z,Y,T).^ 


Note that the factor 11(2:, T, r)/IImax has now the effect of selecting only the resolved 
objects. The number of sources (per solid angle) with a flux larger than S, i.e. the 
cumulative source count distribution N{> S), can be obtained by integrating Eq. (3) 
above S. 

2.2.1. Blazars 

Blazars are among the brightest gamma-ray sources in the sky. They are interpreted 
as Active Galactic Nuclei (AGNs) with the relativistic jet directed towards the observer 
[89, 90]. The emission coming from the jet normally outshines the radiation associated 
with the accretion onto the supermassive Black Hole at the center of the nucleus. The 
Spectral Energy Distribution (SED) is bimodal: the low-energy peak is located between 
ultraviolet (UV) and radio frequencies and it corresponds to the synchrotron emission 
produced by the electrons that are accelerated in the jet. Instead, the high-energy 
peak reaches the gamma-ray band and it is produced by IG of the same population of 
synchrotron-emitting electrons. The seed photons for the IG emission come either from 
the synchrotron radiation itself (i.e., the so-called synchrotron self-Compton) or from 
the accretion disk (external Compton). Blazars can be divided in two categories: BL 
Lacs and Flat Spectrum Radio Quasars (FSRQs). Sources belonging to the first class 
have a nuclear non-thermal emission so strong that the rest-frame equivalent width of the 
strongest optical emission line is narrower than 5 A, while FSRQs have broader lines and 
a spectral index < 0.5 in the radio band. Blazars have also been classified according 
to the frequency of their synchrotron peak, [91, 92]: low-synchrotron-peaked (LSP) 
blazars have a peak in the infrared (IR) or far-IR band, (t's < 10^^ Hz), intermediate- 
synchrotron-peaked (ISP) blazars for i>s in the near-IR or UV (10^^ Hz < us < 10^® Hz) 
and high-synchrotron-peaked (HSP) blazars for us > 10^® Hz, i.e. in the UV band 
or higher. From the study of the 886 blazars present in the Second Fermi LAT AGN 
catalog (2LAC) [93] it was confirmed that FSRQs (which are almost entirely LSPs) are 


®Normally S is the gamma-ray flux above some energy Emin and, therefore, it is related to Feq (z, T, F) 
as S =/p FEo(z,Y,r)dEo. 

-^min 
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generally brighter than BL Lacs. Also, it was possible to determine a correlation between 
the steepness of their energy spectrum (in the gamma-ray band) and the position of the 
synchrotron peak, suggesting that FSRQs (with a low 1 / 3 ) have softer spectra than BL 
Lacs [93]. 

Since the EGRET era, blazars were suspected to play a signihcant role in explaining 
the DGRB emission, with estimates ranging from 20% to 100% [94, 10, 11, 12, 13, 15, 95, 
30, 16]. Fewer AGNs were known before the Fermi LAT (66 in the third EGRET catalog 
[96], mainly FSRQs), preventing reliable population studies to be performed entirely at 
gamma-ray energies. Predictions for the emission of unresolved blazars were obtained 
by means of well-established correlations between the gamma-ray luminosity and the 
luminosity at lower frequencies, either in the radio or in the X-ray band [94, 11]. Thus, 
the strategy followed in these early works was to characterize the blazar population in 
those low-frequency regimes, where the statistical sample was much larger, and then 
export the information up to the gamma-ray range by means of the aforementioned 
correlation between luminosities. 

Ref. [16] considers the correlation between and the luminosity in X-rays, Lx- 
The — Lx connection is a consequence of the relation between the emission of the jet 
(which can be linked to L,y) and the mass accretion onto the supermassive Black Hole 
[97, 98, 99, 100]. This, in turn, correlates with the X-ray luminosity of the accretion 
disk Lx- Ref. [16] considers L.y as the T-parameter in Eqs. (1) and (3). The dependence 
on F in the factor dN/dVdL,ydT is taken care of by the so-called “blazar sequence” 
[101, 102, 103], which predicts how the shape of the blazar SED changes as a function 
of their luminosity. Then, integrating Eq.l over F, the blazar sequence selects, for each 
L-y, the only SED compatible with Ly. On the other hand, dN/dVdLy is the gamma- 
ray luminosity function (LF), ^y{z,Ly). Given the Ly — Lx relation, ^y(z,Ly) can be 
inferred from the X-ray LF 4>x(Lx) ^)- 

<^^{Ly,z) = K^<^xiLx,z)- (4) 

The factor k indicates the fraction of AGNs observed as blazars and a parametrization 
for ^x{z,Lx) is available in Refs. [104, 105, 106]. X-ray AGNs are found to evolve 
positively (i.e. they are more abundant as redshift increases) until a certain redshift peak 
Zc, above which the X-ray LF decreases [104, 105, 106]. R was found that allowing Zc 
to vary with Lx (i.e., what is called a luminosity-dependent density evolution) provides 
a better description of the blazars observed by EGRET than other evolution schemes, 
e.g., pure luminosity evolution or pure density evolution [13]. 

Parameter k in Eq. (4) is left free in the analysis of Ref. [16], together with the 
faint-end slope 71 of the X-ray LF and the proportionality coefficient between Ly and 
Lx-"^ These parameters are determined by a maximum-likelihood fit to the differential 
source count dN/dS of the blazars detected by EGRET. The best-ht model is, then, 
used to determine the contribution of unresolved blazars to the DGRB through Eq. (1). 
Ref. [16] finds that unresolved blazars can explain approximately 45% of the DGRB 
measured by EGRET in Ref. [7] at 100 MeV. 

Similar results are obtained in other works that follow the same formalism: Refs. [13, 


^More precisely, the assumed proportionality is between Lx and the so-called “bolometric luminos¬ 
ity”, i.e. a measure of the total emission power of the source. However, under the hypothesis of the 
blazar sequence, can be uniquely derived from the bolometric luminosity. 
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107, 52] also consider a luminosity-dependent density evolution for the X-ray LF but they 
assume a common power-law energy spectrum at gamma-ray frequencies, instead of the 
SED predicted by the blazar sequence. In Ref. [19] the formalism outlined before is fitted 
to both the blazar dN/dS computed from the IFGL [108] and to the first Fermi LAX 
measurement of the DGRB in Ref. [8] (see also Ref. [109]). The result of the fit proves 
that it is possible to explain both observables with one single population of blazars.® 
Their predictions are re-calibrated in Ref. [68], where the model is fitted against the 
IFGL dN/dS and the Fermi LAT measurement of the DGRB auto-correlation APS 
of anisotropies from Ref. [66] (see Sec. 3). Ref. [68] finds that the DGRB intensity 
energy spectrum and the auto-correlation APS cannot be simultaneously explained in 
terms of unresolved blazars, since a model that fits both the abundance of resolved 
blazars (i.e., their dN/dS) and the DGRB intensity energy spectrum would exceed the 
measured auto-correlation APS. Alternatively, unresolved blazars can reproduce the 
measured anisotropies but, then, their emission only accounts for a maximum of 4.3% 
of the DGRB intensity in Ref. [8] , above 1 GeV. 

Instead of using the — Lx correlation discussed above, Ref. [11] relates to the 
radio luminosity, Lr- The gamma-ray LF is inferred from the radio LF, similarly to what 
done in Eq. (4). The radio LF is taken from Ref. [110]. A power-law energy spectrum is 
assumed at gamma-ray frequencies, while the distribution of spectral indexes dN/dV is 
calibrated to reproduce the sources detected by EGRET. Ref. [11] finds that unresolved 
blazars can ht reasonably well the EGRET DGRB energy spectrum reported in Ref. [6]. 
Other works exploit the relation with similar findings: Ref. [Ill] determines 

the properties of the blazar population by fitting their measured radio LF (assuming 
it follows a pure luminosity evolution), while in Ref. [20] the fit is performed with the 
cumulative source distribution N{> S) of IFGL blazars. The authors of Ref. [20] also 
caution about the use of the sky coverage determined in Ref. [18], since it may be affected 
by systematic uncertainties in the low-flux regime due to low statistics. 

With the advent of the Fermi era, direct population studies at gamma-ray frequen¬ 
cies became possible, without the need to rely on correlations with lower frequencies. 
Following the formalism of Eqs. (1) and (3), in Ref. [18] the Y parameter is taken to 
be the flux S above 100 MeV. The energy spectra of the blazars are assumed to be 
power laws with indexes F characterized by a Gaussian probability distribution inde¬ 
pendent of S. The gamma-ray LF dN/dVdS is combined with the factor dV/dzdQ in 
Eq. (3) into dN/dzdSdLl, which is assumed to be the same for all the sources consid¬ 
ered and to depend on S as a broken power law. No correlation with other frequencies 
is required. The parameters of the broken power law (as well as the mean and stan¬ 
dard deviation of the Gaussian distribution for F) are inferred through a maximum- 
likelihood fit to the blazars in the IFGL with |6| > 20°. The best-fit point has a break 
at 5 = (5.99 ±0.91) x 10“^cm“^s“^ and slopes of 1.58 ±0.08 and 2.44 ±0.11 above and 
below the break [18]. The model is, then, used to determine the cumulative emission of 
unresolved sources. Results indicate that unresolved blazars (with S'min = 0) account 
for 23% of the DGRB detected by Fermi LAT in Ref. [8]. 

Ref. [22] improves the analysis in Ref. [18], attacking one of its main limitations, 
i.e. the fact that the broken power law considered for dN/dzdSdLl is assumed a priori 
and it is not inferred from a specific evolution scheme. In Ref. [22], the analysis is 
restricted to FSRQs: sources are identihed by their (i.e., Y = L^), their redshift and 


®The best-fit point obtained in Ref. [19] favors a 71 much larger than the one found in Ref. [16] for 
EGRET blazars. 


11 




Figure 2: Contribution of unresolved FSRQs to the DGRB as determined by integrating the LF coupled 
to the SED model derived in Ref. [22]. The red hatched band around the best-fit point prediction (solid 
blue line) shows the Ict statistical uncertainty, while the gray band represents the systematic uncertainty. 
Note that, in the figure, the DGRB measured by the Fermi LAT in Ref. [8] (black data points) is refereed 
to as “IGRB”. Taken from Ref. [22]. 


the slope of the power-law fit to their energy spectrum. dN/dVdL^dT is split into the 
gamma-ray LF and the spectral index distribution dN/dV. The former is described by 
a parametric expression inspired by the results in the radio and X-ray bands, while a 
Gaussian distribution is assumed for dN/dT, with no dependence on L.y. A maximum- 
likelihood fit is performed to determine the parameters of the model. 186 FSRQs in the 
IFGL and with |6| > 15° are considered in the fit. Results establish that FSRQs evolve 
positively until a Zc that depends on L^, i.e. a luminosity-dependent density evolution 
performs better than other evolution formalisms. The best-fit model corresponds to an 
unresolved emission which accounts for 9.3('(}g% of the Fermi LAT DGRB in Ref. [8], 
between 0.1 and 100 GeV (see Fig. 2).® The contribution peaks below the GeV scale 
and it becomes more subdominant at higher energies. 

Performing a similar population study for BL Lacs is hampered by the fact that it is 
more difficult to obtain a measure of spectroscopic redshift for these objects due to the 
lack of strong emission lines: indeed, approximately 55% of the BL Lacs in the 2LAC 
do not have an associated z. This issue is somehow alleviated in Ref. [23] by considering 
photometric redshift estimates [113], lower [114] or upper limits on z [113, 114] and 


®In Ref. [22] only sources with L-y > lO^^^erg s“^ are considered when computing the emission of 
unresolved FSRQs. Such a value corresponds approximately to the fainter blazar in the IFGL. Thus, 
the quoted contribution of FSRQs to the DGRB implicitly assumes that no sources are present below 
10 '^‘^erg s“^ and, therefore, it should be considered as a lower limit to the total emission of unresolved 
FSRQs. 
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Figure 3: Diffuse gamma-ray emission from unresolved BL Lacs. Predictions for the best-fit model in 
Ref. [24] are shown embedded in their la uncertainty bands: the summed contribution from LSPs and 
ISPs is plotted by means of the purple dot-dashed line and the purple band, the one from HSP BL Lacs 
by the green band and green dotted line (almost overlapped with the blue band), the one from the sum 
of the two by blue band and solid blue line and the contribution from BL Lacs considered as a unique 
population is depicted by the pink band and dashed pink line. The DGRB data are taken from Ref. [112] 
and are displayed as black points. The gray band and double-dot-dashed black line (orange band and 
dashed red line) represent the cascade emission from the HSP BL Lacs (the whole BL Lac population). 
Taken from Ref. [24]. 


host-galaxy spectral fitting [114]. The 211 BL Lacs studied in Ref. [23] are taken from 
Ref. [18] and analyzed by means of the same pipeline applied in Ref. [22] to FSRQs: the Y 
parameter in Eq. (1) is again but the mean of the Gaussian distribution of spectral 
indexes depends now linearly on log]^Q(L-y). Only sources with > 7 x lO^^erg s“^ 
are considered. The best-fit model suggests that the number density of faint BL Lacs 
(probably HSPs) decreases with redshift, while BL Lacs with > 10^^’®erg s“^ are 
more numerous at large redshift, i.e. more similar to the positive evolution FSRQs. 
The redshift estimates adopted for the BL Lacs without spectroscopic information are 
crucial to constrain the evolution of their LF. Results indicate that unresolved BL Lacs 
contribute to 7.71^3% of the Fermi LAT DGRB in Ref. [8], between 0.1 and 100 GeV. 
Due to the large density of low-luminosity hard sources at low redshift, the emission of 
unresolved BL Lacs is expected to be harder than that of FSRQs and, thus, BL Lacs 
may play a more significant role at higher energies. 

This hypothesis was tested by Ref. [24], where the authors considered a set of 148 
BL Lacs with redshift and synchrotron peak frequency obtained from the 2FGL 
[115]. Their model for the gamma-ray LF is fitted to the observed cumulative source 
distribution A^(> S), the gamma-ray LF and the redshift distribution dN/dz of the 
detected sources. The energy spectra of the sources are obtained from a combination 
of Fermi LAT and Imaging Air Gherenkov Telescopes data. A luminosity-dependent 
density evolution is found to provide the best ht to the data. Pure power laws, log- 
parabolae and power laws with exponential cut-offs are considered as possible SEDs, 
with the last one corresponding to the most accurate description of the BL Lacs in the 
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Figure 4: Diffuse emission arising from blazars (with or without EBL absorption), in comparison with 
the intensity of the total emission from sources (both resolved and nnresolved), called here “EGB” (red 
data points, from Ref. [9]). Taken from Ref. [25] 


sample. The sources were considered as either one single population, or split into HSPs 
and a second sub-class including ISPs and LSPs. In their best-fit model, HSPs dominates 
the dN/dS below S' = 5 x 10“®cm“^s“^ and their SED extends to much higher energies 
than in the ISP-I-LSP class (the best-fit cut-off energy is 910 GeV for HSPs and 37 GeV 
for the class of ISPs and LSPs). That is the reason why the cumulative emission from 
HSPs (computed from Eq. (1) above L-y > lO^^erg s“^) can extend up to very high 
energies and it is able to explain the whole DGRB emission reported in Ref. [112] above 
few tens of GeV (see Fig. 3). Between 0.1 and 100 GeV, unresolved BL Lacs account 
for ~ 11% of the Fermi LAT DGRB in Ref. [112], in agreement with Ref. [23]. 

Ref. [25] repeated the analysis of Ref. [23] on a sample of 403 blazars from IFGL, 
this time considering both FSRQs and BL Lacs as one single population by allowing 
the spectral index distribution to depend on A double power-law energy spectrum, 
proportional to [(Eo/Ef,)^-’^ -|- (Eo/Eb)^'®]“^, is assumed and the energy scale E;, is found 
to correlate with the index P obtained when the SED is fitted by a single power law. 
The same LF used in Ref. [23] and based on a luminosity-dependent density evolution 
is implemented in Ref. [25], together with other evolution schemes. They all provide an 
acceptable description of the blazar population, even if the luminosity-dependent density 
evolution is the one corresponding to the largest log-likelihood. The predicted cumula¬ 
tive emission of blazars (FSRQs and BL Lacs, resolved and unresolved) can be seen in 
the Fig. 4 as a dotted blue band, compared to the total emission from resolved and unre¬ 
solved sources taken from Ref. [9] (labeled “EGB” here, red data points). Blazars (both 
resolved and unresolved) accounts for the 50(](}^% of the total emission from resolved 
and unresolved sources, above 100 MeV. Unresolved blazars, on the other hand, are 
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responsible for approximately 20% of the new Fermi LAT measurement of the DGRB 
in Ref. [9], in the 0.1-100 GeV energy band. 

2.2.2. Misaligned Active Galactic Nuclei 

Under the AGN unification scenario, the parameter that discriminates among dif¬ 
ferent classes of AGNs is the viewing angle [89]. A value of 14° separate blazars (with 
the jet pointing towards the observer) from non-blazars, i.e. MAGNs. Among MAGNs 
it is possible to further distinguish between radio galaxies (with a viewing angle larger 
than 44°) and radio quasars [116]. Radio galaxies are classihed either as Fanaroff-Riley 
Type I or Type II (FRI and FRII, respectively) according to their morphology [117]. 
The emission of FRIs peaks at the center of the AGN and it is dominated by two-sided 
decelerating jets. On the other hand, FRIIs are brighter and they are characterized by 
edge-brightened radio lobes with bright hotspots, while jets and core (when detected) are 
subdominant. FRIs and FRIIs are normally interpreted as the misaligned counterparts 
of BL Lacs and FSRQs, respectively. Indeed, the detection of MAGNs in the IFGL 
and 2FGL conhrms that, similar to the case of BL Lacs and FSRQs, FRIs have harder 
spectra than FRIIs [118]. 

The mechanisms of gamma-ray production in MAGNs are less clear than for blazars: 
the SED exhibits a bimodal structure, with the low-frequency peak due to synchrotron 
radiation, while at higher energies, the bulk of gamma-ray emission may be due to syn¬ 
chrotron self-Compton [119]. External IC may also contribute [120]. If present, external 
IC would be localized within one parsec from the center of the AGN, while synchrotron 
self-Compton emission is produced outside that region [119, 120]. Furthermore, spa¬ 
tially coincident with the radio lobes, there may also be emission from IC off the Cosmic 
Microwave Background (CMB) [27, 26] (see the discussion at the end of this section). 

With no Doppler boost, MAGNs are expected to be less bright but more abundant 
than blazars, making them potentially important contributors to the DGRB [28]. 15 
MAGNs were reported in Ref. [118], precluding the possibility of deriving their gamma- 
ray LF directly from gamma-ray observations, as done with blazars. Yet, the LF of 
MAGNs is well known at radio frequencies. As done in the previous section, once a 
correlation between their radio and gamma-ray luminosities is established, it is possible 
to deduce the gamma-ray properties of MAGNs by studying the sources in the radio 
band. 

Ref. [28] considers L.y (between 0.1 and 10 GeV) for 10 MAGNs detected by the Fermi 
LAT in Ref. [118]. The author studies the possibility of a linear correlation between 
log;^o(-^ 7 ) logio(-^r’)! where is the radio luminosity. The slope of the observed 
correlation is very similar to the one found by Ref. [121] for blazars. Then, $..y(z,L^) 
is determined in terms of the radio LF, <hr( 2 ;,Lr), as <h^(z,L^) = K,^r{z,Lr)dLr/dL.y. 
This is similar to Eq. (4) for blazars. Ref. [28] takes 4>r from Ref. [122] and the value of 
K is tuned to reproduce the number of sources observed by the Fermi LAT. All MAGNs 
are assumed to share the same energy spectrum, i.e. a power law with an index of 2.39. 
Such a value is the mean spectral index among the 10 MAGNs considered. The emission 
from unresolved MAGNs is, finally, estimated from Eq. (1), for L.y > 10^®erg s The 
result indicates that ~ 25% of the DGRB measured by Fermi LAT in Ref. [8] above 100 
MeV can be explained in terms of unresolved MAGNs. 

Ref. [29] improves the analysis in Ref. [28] by properly estimating the uncertainties 
involved: the authors consider L.y (defined above 100 MeV) of 12 MAGNs in the 2FGL 
and the correlation with at 5 GHz. They make use of the total radio luminosity, 
Lr,tot, and the so-called “radio core luminosity” Lr,core, defined as the emission from 


15 


the central arcsecond-scale region of the source. Both the log;^o(-^ 7 ) “ logig(Lr,core) and 
log;^o(-^ 7 ) “ relations are considered. The gamma-ray LF is inferred from 

the radio LF as follows: 


*^7(2:, L..y) — 


logio(L^,0 

loglo(-^^7) ’ 


(5) 


where i stands for “total” or “core”.^*^ If Lr, tot is used, <hr,tot is taken from Ref. [122]. On 
the other hand, it is not possible to build the core radio LF directly from radio observa¬ 
tion, due to the scarce data available [123]. In Ref. [29], a relation is assumed between 
Lr,core and Lr^tot [124], so that d'r.core Can be derived from The factor Ki in Eq. (5) 

is tuned to reproduce the number of observed MAGNs. Different values for k are ob¬ 
tained in Ref. [29], depending on which LF is considered (total radio or core radio) and 
on how the uncertainties in the log^g(Lg,) — log]^Q(Lr_j) correlations are treated. Finally, 
the emission from unresolved MAGNs is computed as in Eq. (1), assuming a Gaussian 
distribution for the spectral indexes. The results (for > lO^^erg s“^) are summarized 
in Eig. 5. For the best-fit model, the contribution of MAGNs accounts for ~ 25% of 
the Fermi LAT DGRB in Ref. [8] above 100 MeV. The value agrees very well with that 
given in Ref. [28]. Yet, in Ref. [29], the prediction is embedded in an uncertainty band 
with a size of almost one order of magnitude. Such large uncertainty mainly comes from 
the possible different values for n. 

Given the complex morphology of the emission in radio galaxies, it is possible that 
different regions in the source emit gamma rays which are not accounted for by the anal¬ 
yses presented above. Ref. [27] considers a scenario in which the non-thermal electrons 
responsible for the synchrotron emission in the radio lobes of FRIIs could emit gamma 
rays by IC off the CMB. Indeed, gamma-ray emission spatially coincident with the ra¬ 
dio lobes of Gentaurus A has been recently observed by the Fermi LAT [125]. Such a 
component would peak at around 1 MeV and could explain up to 10% of the DGRB in 
Ref. [8] below 1 GeV. 

Similarly, Ref. [26] considers the X-ray-emitting kpc-scale jets of FRIs as possible 
contributors to the DGRB. The established synchrotron origin of those X-rays suggests 
that the same non-thermal electrons could also emit gamma rays through IC with the 
ambient photons of the host galaxy. However, even assuming this is common for all 
FRIs, the contribution adds up to only ~ 1% of the DGRB reported in Ref. [7]. 


2.2.3. Star-forming galaxies 

The majority of the gamma-ray emission detected by the Fermi LAT is associated 
with the MW [79]. This diffuse Galactic foreground is produced by the interaction of 
Galactic CRs (mainly protons and electrons) with the Galactic interstellar medium and 
interstellar radiation field. Other SFGs similar to the MW are expected to shine in 
gamma rays thanks to the same emission mechanisms. Up to now, the Fermi LAT has 
only detected a few SFGs other than the MW: M31 and M33 [126], the Large and Small 
Magellanic Clouds [127, 128] in the Local Group and the Circinus Galaxy [129], M82, 
NGC 253 [130, 131], NGC 1068 and NGC 4945 [33]. Ref. [132] also reported the detection 
of gamma-ray emission from NGC 2146. Being intrinsically faint but numerous, SFGs 
are expected to contribute significantly to the DGRB. 

Massive stars in SFGs emit the majority of their light in the UV band. The emission 
is then absorbed by interstellar dust and re-emitted as IR light. The IR luminosity can be 


^°The luminosity functions in Ref. [29] are defined per units of log(L), hence the factor 
logiofLr.ij/logiofL^) in Eq. (5) as compared to Eq. (4). 
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Figure 5: Diffuse gamma-ray flux due to MAGNs as a function of the gamma-ray energy. The dashed 
black line is obtained from the best-fit to the Lr.core — Lj correlation. The cyan shaded area is derived 
by allowing for the uncertainty on the Tr.core — correlation and by including all the models that are 
less than Icr away from the best-fit point. The green band corresponds to the case with ac = 1 at Icr 
confidence level. The red dot-dashed cnrve shows the diffnse flux obtained when assuming a correlation 
between the total radio luminosity and the gamma-ray luminosity. The black squares corresponds to the 
DGRB measured by the Fermi LAT in Ref. [8] and the magenta dashed curve represents the best-fit to 
the data. Taken from Ref. [29]. 


used as a good tracer of the star-forming rate (SFR) t/;, i.e. the amount of mass converted 
in stars per unit time [133]. The same massive stars finally explode into core-collapse 
supernovae, leaving behind supernova remnants which are considered to be the main 
sources of accelerated CRs on galactic scales [134], The leptonic component of these CRs 
is responsible for the synchrotron radio emission observed from SFGs and they may also 
contribute at gamma-ray frequencies when interacting with the interstellar radiation held 
of the host galaxy (through IC or bremsstrahlung). However, such leptonic gamma-ray 
emission is expected to be subdominant with respect to the hadronic one [35] , since CR 
protons have an intrinsically larger injection rate than CR electrons [135, 34]. Hadronic 
gamma-ray emission mainly comes from inelastic interactions of CR protons with the 
nuclei of the interstellar medium, producing neutral pions that decay into gamma rays 
with an energy spectrum that peaks at around 300-400 MeV. 

The so-called “initial mass function” determines the relative number of stars pro¬ 
duced in any star-formation event in a galaxy, as a function of the stellar mass [136, 137]. 
Under the assumption of an universal high-mass end of the initial mass function, the 
SFR of a galaxy is proportional to the rate of supernova explosions and, thus, to the 
CR abundance and the CR-induced emission [138, 34], Then, it is expected that, in a 
generic galaxy, the gamma-ray and radio luminosities (both associated with CRs) are 
correlated with the IR emission, which depends on the SFR. 
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The cumulative gamma-ray flux produced by IC of CR electrons in all SFGs at all 
redshifts is computed in Ref. [35]. In this work, the emission of a generic galaxy is 
modeled based on a template which is tuned to reproduce the emission of the MW. The 
cumulative gamma-ray flux depends on the abundance of SFGs. Ref. [35] assumes that, 
at a certain redshift, the total emission scales as the so-called cosmic SFR Pi,{z), i.e. the 
mass converted to stars per unit time and comoving volume: 

( 6 ) 

where dN/dVd'ip is the comoving density of galaxies per unit SFR. Parametric fits for 
Pi,{z) are available from the observation of the total luminosity density at various wave¬ 
lengths [139, 140, 141]. Ref. [35] finds that the IC-induced emission is always subdom¬ 
inant with respect to the hadronic one. The two become comparable only above 100 
GeV. 

Indeed, the majority of the studies on the SFGs focus on their hadronic emission 
and they neglect the contribution from primary electrons. The hadronic gamma-ray 
luminosity of a SFG, as a function of the observed energy Eq, can be written as follows: 


L-yiEo) = j Ur dV = F.,r0^.77(-^em) (7) 

where Flem is the energy in the emitter frame, is the pionic gamma-ray produc¬ 

tion rate per interstellar hydrogen atom and nn is density of hydrogen atoms in the 
interstellar medium [31, 20]. Integrated over the volume of the medium, A^h gives the 
total number of hydrogen atoms, which can be expressed in terms of the total interstellar 
gas mass Mgas in the galaxy as XHMgas/uip. Xr ~ 0.7 is the hydrogen mass fraction 
and rup is the mass of the proton. As for in Eq. (7), it is the product of the flux 

of GR protons (averaged over the volume of the galaxy) and the cross section for the 
production of gamma rays [142, 143]. 

The distribution and propagation of CR protons is governed by the diffuse-loss equa¬ 
tion, which depends on their injection rate, diffusion, energy losses and possible inelastic 
interactions (see Ref. [34] for a recent review). Different scenarios are possible depending 
on the strength of those terms. Here we only mention two possibilities that provide good 
descriptions to different typologies of galaxies. The first is the so-called escape regime 
and it corresponds to a situation in which the energy losses of CR protons are dominated 
by their escape from the diffuse region of the galaxy. An equilibrium is reached between 
the energy losses and the acceleration of CR protons so that their flux is proportional to 
the product of the SFR of the galaxy and the CR escape path-length, Aesc: ‘bp oc AescV’- 
Consequently, oc [144, 145, 31]. The escape regime provides a good descrip¬ 

tion of the diffuse foreground of the MW. This motivates the use of the gamma-ray 
production rate of the MW, to normalize the proportionality relation between 

and Then Eq. (7) becomes: 


L-y ipj) , Afgas) 


ArF 


MW 

77^^77 


-^gas 

m-p V'Mw’ 


( 8 ) 


where V'mw is the SFR of the MW. This relation provides a good description of the 
so-called quiescent or normal SFGs, characterized by properties similar to those of the 
MW. 

The second scenario considered for the modeling of SFGs is the so-called calorimetric 
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regime: in this case the energy losses of CR protons is mainly due to their inelastic 
interactions. This means that protons lose all their energy into pions and SFGs act 
effectively as calorimeters. Their gamma-ray lumonisity, then, can be computed from 
the amount of energy available to CRs. Under the paradigm that supernova remnants 
are the primary source of CRs, the energy available in the form of CRs for a calorimetric 
SFC is proportional to the supernova rate, multiplied by the energy FigN released per 
supernova and by the fraction ry of that energy going into CR protons (i.e. the so-called 
acceleration efficiency) [138, 33]. In turn, the supernova rate is proportional to the SFR, 
so that, finally, 

L^{ip) (xip Esni]. (9) 

Starburst galaxies are well modeled as proton calorimeters: normally brighter and less 
numerous than quiescent galaxies, starburst ones are characterized by at least one region 
undergoing intense star formation. This is often induced by major merger events or by 
bar instabilities, leading to a large gas density [146]. 

Several works analyze quiescent and starburst SFCs separately. Refs. [144, 145, 31], 
e.g., focus on quiescent galaxies. Their contribution to the DCRB is computed by using 
a formalism similar to Eq. (1), where is considered as the T-parameter. Then, by 
means of Eq. (8), is written as a function of Mgas and i/j. In Refs. [144, 145], the 
former is factorized out of the integral in dL^ in Eq. (1), assuming that Mgas only 
depends on redshift. The integration in dL^ now only depends on the SER and it 
amounts to the cosmic SFR p-k{z) in Eq. (6). The net result is that the cumulative 
emission of quiescent SFCs scales with redshift as the product of Pi,{z) and the average 
amount of gas. Refs. [144, 145] take Pi,{z) from Refs. [139, 140], respectively, and derive 
the average gas mass assuming that the total baryonic mass of a galaxy (stars and gas) 
remains constant time and that it can be computed de-evolving backwards the cosmic 
history of the SFR density [144, 145]. They find that quiescent SFCs can contribute 
significantly to the DGRB, especially below 1 GeV, but they are only subdominant at 
higher energies. 

Ref. [31] assumes the same formalism than Ref. [144] but it relates the gas content 
of a galaxy to its SFR, using the following relation: 

Mgas(iA, ^) = 2.8 X lO^Mo (1 + z)-^-^'^^ ( ^Jr-I ) ' 

This is obtained by the so-called Kennicutt-Schmidt law [147, 148], which links the 
surface density of star formation to that of the gas. Then, the integration over in 
Eq. (1) is converted into one in SFR. In turn, the SFR is related to the IR luminosity, 
Lir, assuming a direct proportionality between the two [149]. The contribution of SFGs 
now depends on the IR LF, <hiR. Ref. [31] considers the IR LF from Ref. [150]. Results 
are provided for both a pure luminosity evolution and a pure density evolution for 
<I>IR. The difference between the two prescriptions amounts to approximately a factor 
of 2 (see the dashed and long-dashed lines in Fig. 6). An intermediate scenario that 
mediates between the two evolution schemes indicates that unresolved SFGs account for 
~ 50% of the DGRB intensity reported in Ref. [8] below 10 GeV. Above that energy, 
the contribution of SFGs goes down rapidly due to the softness of the energy spectrum 
assumed, i.e. a power law with a slope of 2.75. 

Ref. [20] proposes two alternatives to estimate Mgas: in the first one, the gas mass 
is not obtained from Eq. (10) but taken to be some fraction of the stellar mass M*. 
The latter can be constrained by combining the findings of Refs. [151, 152]. In the 
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second scenario, the mean gas mass is assumed to scale with the cosmic SFR divided by 
a parameter linking the SFR to the density of clouds of molecular hydrogen [153, 20]. 
The two methods differ by approximately a factor of 10 in their predictions for the 
contribution of SFGs to the DGRB. In the second one, SFGs are able to account for the 
whole DGRB of Ref. [8] between 200 MeV and 1 GeV. 

Starburst SFGs have also been considered: Ref. [138] works under the hypothesis that 
these objects are in the calorimetric regime. From Eq. (9), their gamma-ray luminosity 
is proportional to the SFR which, in Ref. [138], is related to the IR luminosity. Thus, 
the total gamma-ray emission from unresolved starburst SFGs can be obtained simply 
by rescaling the total diffuse extragalactic IR background taken, e.g., from Ref. [154]. 
Unresolved starburst galaxies are found to account for ~ 10% of the DGRB detected 
by EGRET in Ref. [7], at least at the GeV scale (see the double-dotted-dashed line in 
Fig. 6). The main uncertainty affecting this prediction stems from the unknown star- 
burst fraction, i.e. the fraction of the IR background emission associated with starburst 
galaxies in the calorimetric regime. Observations suggests that such a fraction is quite 
low at 2 ; = 0 but that it can approximate to I at earlier epochs, when the star formation 
is much more efficient [155, 156]. However, see also Refs. [157, 20, 158, 34] for different 
values. 

Both quiescent and starburst SFGs are modeled at the same time in Ref. [32]. The 
relations between Lj and ipMgas and between and ip for quiescent and starburst 
SFGs respectively (see Eqs. (8) and (9)) are determined by fitting the results of 4 SFGs 
detected by the Fermi LAT. On the other hand, the abundance of galaxies is inferred 
from simulated galaxy catalogs [159, 160]. Starburst SFGs are found to be always 
subdominant and the total (quiescent and starburst) gamma-ray emission is between 
5.4% and 9.6% of the DGRB intensity reported by the Fermi LAT in Ref. [8] above 100 
MeV (see the black solid line in Fig. 6). 

In Ref. [34] the authors model the emission from the MW and from M82, which are 
considered as templates of quiescent and starburst SFGs, respectively. The results are 
used to determine for the two typologies of SFGs. Their total emission is assumed 
to evolve with redshift following p-k{z), modulated by functions fi{z) (where i stands for 
either “starburst” or “quiescent”) that fix the relative abundance of the two sub-classes. 
They find that SFGs can be responsible from 4% to 76% of the DGRB of Ref. [8] in 
the GeV range and that their contribution cannot reproduce the data below the GeV 
or above few tens GeV. In their fiducial model quiescent SFGs always dominate over 
starbursts. 

Ref. [33] follows an alternative approach to compute the emission of unresolved SFGs. 
The authors consider the 8 SFGs detected by the Fermi LAT and 64 galaxies observed 
in radio and IR but for which only upper limits are available in the gamma-ray energy 
range. These are used to determine a correlation between log;^o(^ 7 ) (defined between 
0.1 and 100 GeV) and the IR luminosity logig(LiR,) (defined between 8 and 1000 pm). 
Similarly to what done for blazars and MAGNs, this correlation is used to infer the 
properties of SFGs in the gamma-ray band from the study performed at IR frequencies. 
In particular, it is assumed that the gamma-ray LF can be written in terms of <biR, as 
= 4>iRdlogio(LiR)/dlogio(L7). 

The contribution to the DGRB, then, is computed following Eq. (1), with Y = 

The IR LF is measured in Ref. [162] from data gathered by the Spitzer Space Telescope. 


done for MAGNs, gamma-ray and IR LF are defined here per unit of log{Lj) and of log(LiR), 
respectively. 
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Figure 6: Estimated contribution of unresolved SFGs (both quiescent and starburst) to the DGRB 
emission measured by the Fermi LAT (black points, taken from Ref. [8]). Two different spectral models 
are used: a power law with a photon index of 2.2 (red line), and a spectral shape based on a numer¬ 
ical model of the global gamma-ray emission of the MW [135] (white line, inside gray band). The 
shaded regions indicate combined statistical and systematic uncertainties. Several other estimates for 
the intensity of unresolved SFGs are shown for comparison. For the double-dotted-dashed line starburst 
galaxies are treated as calorimeters of CR nuclei as in Ref. [138]. Ref. [31] (dashed and long-dashed 
line) considers the extreme cases of either pure luminosity and pure density evolutions. The solid black 
line shows predictions from Ref. [32], obtained from simulated galaxies and semi-anal 3 rtical models of 
galaxy formation. Two recent predictions from Ref. [20] are plotted by dotted and dash-dotted black 
lines: the former assumes a scaling relation between IR and gamma-ray luminosities and the latter uses 
a redshift-evolving relation between the gas mass and the stellar mass of a galaxy. Taken from Ref. [33]. 


In Ref. [33] only objects below z = 2.5 are considered, since <I>ir is not well determined 
for higher redshifts. The Lin ~ -^7 correlation is assumed to stand valid up to that 
redshift, even if only SFGs below z ~ 0.05 are used in Ref. [33] to derive it. Starburst 
galaxies are observed to have a harder energy spectrum than quiescent ones [163, 130], 
at least in the Local Group. Assuming that this is a general property, the steepness 
of the contribution of unresolved SFGs to the DGRB will depend on the fraction of 
starburst galaxies, which is unknown. Ref. [33] assumes that all SFGs share the same 
energy spectrum and it considers two extreme cases: a power-law with a slope of 2.2 
(typical of starburst galaxies) and the spectrum of the diffuse foreground of the MW 
[135] , which reproduces well the behavior of quiescent SFGs. Results are summarized in 
Fig. 6 by the red and gray bands, respectively. The figures show that unresolved SFGs 
can explain between 4% and 23% of the DGRB measured by the Fermi LAT in Ref. [ 8 ] 
above 100 MeV. 

The results of Ref. [33] have been updated in Ref. [161], by improving the modeling 
of Thanks to the recent detection of a larger number of high-^; sources by the 

Herschel Space Observatory [164, 165, 166, 167], it is now possible to extend the LF 


21 




















Figure 7: Diffuse gamma-ray emission as a function of the energy (without EBL attenuation) for the 
three different SFG contributions to the DGRB: quiescent SFGs (green, labeled “NG” in the figure), 
starburst SFGs (blue) and SFGs hosting an AGN (black). The solid pink indicates their sum. The 
DGRB from Ref. [8] is plotted in red. Taken from Ref. [161]. 


up to 2 ; ~ 4 and to consider separately the contribution of different classes of SFGs. In 
particular, Herschel classifies its sources into quiescent galaxies, starburst ones and SFGs 
hosting an obscured or low-luminosity AGN.^^ At z = 0, starburst SFGs are found to be 
subdominant with respect to the other two classes, but the three families have different 
evolution and, by z = 1, the number of AGN-hosting SFGs is approximately twice the 
number of quiescent and starburst galaxies combined. Ref. [164] provides analytical 
fits to the LF of the three classes. From these, the emission of unresolved SFGs can 
be computed, assuming a broken power law as energy spectrum. The results are 
summarized in Fig. 7 and they show that unresolved SFGs are responsible for ~ 50% of 
the Fermi LAT DGRB of Ref. [8] in the range between 0.3 and 30 GeV. Note that the 
emission is dominated by SFGs hosting AGNs. 

2.2.4- Millisecond pulsars 

Pulsars are highly magnetized and rapidly spinning neutron stars, with a beam of 
radiation that periodically intersects the Earth. Their initial spin P decreases due to 


note that, even if these SFGs host an AGN, their IR luminosities are expected to be dominated 
by the star-forming activity rather than the one associated with the AGN and, therefore, Lir remains 
a useful tracer of the SFR. 

^®The low-energy slope is fixed to -1.5 in all cases, while the high-energy one is assumed to be -2.7 
for quiescent galaxies and -2.2 for starburst ones. The class of AGN-hosting objects is additionally split 
into “starburst-like” and “quiescent-like” galaxies, with slopes of -2.2 and -2.7, respectively. 
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magnetic dipole braking [168], so that the time derivative of the period P can be written 
as follows [36, 37, 169, 170, 38]: 

F = 9.8xl0-(§)'(f) (11) 

where B is the surface magnetic field. The loss of kinetic energy associated with the 
slowing down of the spinning, i.e. the so-called “spin-down luminosity” E, is 

E = i7r^M^, (12) 

where M is the moment of inertia of the neutron star. The spin-down luminosity is con¬ 
verted, with some efficiency, into radiation. MSPs were traditionally detected in radio, 
while, nowadays, thanks to the Fermi LAT, an increasingly large number of objects is 
observed also in the gamma-ray band. The shape of their gamma-ray energy spectra 
suggests that the emission comes from curvature radiation. This is a mechanism similar 
to synchrotron radiation, in which gamma rays are produced by relativistic charged par¬ 
ticles following the curved force lines of a magnetic field [171, 172, 173]. Self-synchrotron 
Compton possibly also contributes [174]. 

Pulsars are classified in terms of their period and sources with P < 15 ms are 
referred to as MSPs. These are generally part of a binary system and their higher spin 
is the result of the large angular momentum transferred from the companion object 
[175, 176, 177, 178, 179]. Thanks to lower surface magnetic fields, MSPs have smaller 
P and are, therefore, older than “normal” (i.e., young) pulsars. A longer life cycle may 
compensate the intrinsic lower birthrate so that MSPs are expected to be as abundant 
as normal pulsars [180]. Moreover, having had the time to orbit many times around 
the Galaxy, the distribution of MSPs is expected to be uncorrelated with their birth 
locations, potentially extending to high Galactic latitudes [36]. 

The conversion of energy loss E into is parametrized by an empirical relation of 
the form 

= r]E° (13) 

[181], where r] is the called “conversion efficiency”. The case of a = 0.5 is motivated in 
Refs. [182, 183, 181] and adopted in Refs. [184, 37, 169], even if some of those references 
also consider the case of a = 1. Assuming that Eq. (13) is valid both for young pulsars 
and for MSPs, their luminosities are related by 

^ -Pmsp f Pmsp \ 
ryoung Ip I ) 

-'young \-'Ryoung/ 

for a = 1. Typical values are Pmsp = 3 ms, Pyoung = 0.5 ms, i\isp = 10“^® and 
Ryoung = 10“^®. These suggest that MSPs can be brighter than normal pulsars, by 
a factor of few tens. It is, therefore, reasonable to speculate that these sources can 
contribute significantly to the DGRB emission [36, 37, 185, 169, 38]. On the other 
hand, the cumulative emission of normal pulsars would be quite anisotropic and mainly 
localized along the Galactic plane, thus hardly compatible with the DGRB [186, 184]. 

The second Fermi LAT catalog of gamma-ray pulsars (2FPG) contains 117 sources, 
40 of which are MSPs [170]. The number is too low to build a MSP LF only from Fermi 
LAT data. Also, given the uncertainties on the mechanisms of gamma-ray emission, it is 
not possible to postulate correlations with other frequencies, as done for blazars, MAGNs 
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Unresolved MSPs flux In the high-latitude region 



E [GeV] 

Figure 8: Prediction for the contribution of unresolved MSPs to the DGRB, as derived from 1000 Monte 
Carlo simulations of the MSP population of the MW. The red solid line represents the mean spectrum 
distribution (see Ref. [38] for further details), while the orange band corresponds to its la uncertainty 
band. The black points refer to the preliminary DGRB measurement anticipated in Ref. [112]. Taken 
from Ref. [38]. 


and SFGs. This is the reason why Refs. [36, 37, 169, 38] rely on Monte Carlo simulations 
in order to properly describe the population of MSPs. Probability distributions of the 
most relevant quantities (e.g., radial and vertical position of the source, its period and 
surface magnetic field) are derived from the pulsars detected in radio [187, 180, 188]. At 
date, the largest catalog of pulsars is the Australia Telescope National Facility Pulsar 
catalog, containing 1509 sources, out of which 132 are MSPs [189]. Ref. [38] analyses 
the objects in the catalog and establishes that they are well described by the following 
prescriptions: i) a Gaussian distribution for the radial distance, R, from the center 
of the MW projected on the Galactic plane, i.e. dN/dR oc exp[—(i? — {R))/Rq], ii) 
an exponential distribution for the vertical distance, z, from the Galactic plane, i.e. 
dN/dz oc exp(— z/zq) [187, 184], in) a log-Gaussian distribution for the period P (in 
contrast to previous works [36, 37, 169]) and iv) a log-Gaussian distribution for the 
magnetic field. 

Ref. [38] also studies the relation between and E. The data, however, are affected 
by quite large errors which prevent a statistically meaningful fit to the data. Thus, 
Ref. [38] identihes a benchmark case with a = 1 that provides a qualitative good de¬ 
scription of the data, and an uncertainty band that encompasses reasonably well the 
distribution of the data points. This completes the characterization of the MSP popu¬ 
lation and synthetic sources can be generated by randomly drawing from the assumed 
distributions. Each source is labeled as “resolved” or “unresolved” depending on whether 
its flux is larger or smaller than the sensitivity of the telescope at the position of the 
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source in the sky. Simulated sources are accumulated until the number of “resolved” 
ones equals the amount of detected MSPs: Refs. [184, 37] consider the MSPs detected 
in radio in the Australia Telescope National Facility catalog [190], while Refs. [169, 38] 
the sources detected by the Fermi LAT.^^ This calibrates the Monte Carlo data and 
it allows the computation of the MSP contribution to the DGRB simply by summing 
over the “unresolved” sources. Their energy spectrum is hxed to a power law with an 
exponential cut-off, i.e a functional shape that reproduces well the MSPs in the 2FPC. 
The slopes and cut-off energies are assumed to have Gaussian distributions. The results 
of Ref. [38] are shown in Fig. 8, from which it is evident that MSPs are only a subdom¬ 
inant component to the DGRB, responsible for less than 1% of the intensity measured 
by the Fermi LAT in Ref. [8]. This value is lower than the previous estimates from 
Refs. [37, 169], in which a different modeling of the MSP population was adopted. 

2.2.5. Other astrophysical components 

Several other potential contributors to the DGRB have been proposed in the past. 
Some correspond to unresolved astrophysical populations not considered in the previous 
sections, while others are intrinsically diffuse processes. In this section, we focus on 
astrophysical scenarios, while a potential DM-induced emission will be discussed in detail 
in Sec. 2.3. 

• Clusters of galaxies: it is believed that huge amounts of energy, of the order of 
lO^i _ 10®^ erg, is dissipated in the shocks associated with the assembly of clusters 
of galaxies [191, 192]. A fraction of this energy can go into the acceleration of 
GRs, even if the details of the acceleration mechanisms are still uncertain [193]. 
Accelerated GRs would produce gamma rays by means of i) the decay of pions 
produced by the interaction of GR protons with the intracluster medium and ii) 
IG scattering of primary GR electrons or of the electrons produced by the pion 
decays in point i). Yet, no gamma-ray emission has been detected from galaxy 
clusters. On the other hand, the interpretation of the observed radio emission as 
synchrotron radiation confirmed the presence of accelerated electrons and magnetic 
fields [194, 195]. However, not all known clusters emit in radio [196, 197] and it is 
not clear why some objects are radio-quiet. 

The IG-induced gamma-ray emission [198, 199] depends on the abundance of GR 
electrons, on the intensity of the magnetic held and on the so-called acceleration 
efficiency i-e. the fraction of the thermal energy density produced by the shocks 
that is transferred to GR electrons. A Halo Model [200, 201] , linking the properties 
of the shock to those of the DM halo hosting the cluster [202, 203], is often used 
to predict the IG-induced emission of clusters. Results can be tested with N-body 
simulations. An acceleration efficiency of = 0-05 is typically assumed, following 
similar values measured in shocks surrounding supernova remnants. For this 
Ref. [202] hnds that the cumulative IG-induced gamma-ray hux from unresolved 
clusters can explain, at most, 10% of the EGRET DGRB in Ref. [6]. Similar 
values are obtained by Refs. [204, 205, 206, 107, 207]. Yet, the non-detection of 
gamma rays from the observation of the Goma galaxy cluster suggests even lower 
efficiencies (^e < 1%) [208], which would further decrease the contribution of this 
source class to the DGRB. 


this case the sensitivity of the Fermi LAT is taken from Ref. [170]. 
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On the other hand, the gamma-ray emission expected from pion decay can be 
estimated as a function of the amount of gas in the intracluster medium, of the 
injected spectrum of CR protons and of the intensity of the magnetic fields [209]. 
Ref. [210] estimates that this hadronic gamma-ray emission can only account for 
less than few percents of the DGRB reported by EGRET in Ref. [6]. This result 
was later reduced to less than 1% of the Fermi LAT DGRB in Ref. [8], after that 
mock catalogs of clusters from Ref. [211] were considered in Ref. [208]. 

The most recent predictions for the emission of unresolved clusters come from 
Ref. [39]. In this paper, a correlation between the mass of the cluster and its 
gamma-ray luminosity is assumed and it is calibrated to reproduce the number 
of clusters detected in radio during the Radio Astronomy Observatory Very Large 
Array sky survey [212]. It is also required that results are compatible with the non¬ 
detection of the Coma cluster by the Fermi LAT [208] and of the Perseus cluster 
by MAGIC [213]. A second, more physically motivated, model for the distribution 
of CRs and of the intracluster medium distribution is also considered in Ref. [39]. 
In this case, the intracluster medium is reconstructed from X-ray observations 
and the CR spatial and spectral distributions are based on hydrodynamic iV-body 
simulations [211]. The two scenarios agree in finding that the contribution of 
galaxy clusters to the DGRB is negligible. 

• Interactions of ultra-high-energy cosmic rays with background radia¬ 
tion: UHECRs, with energies larger than 6 x 10^® eV are attenuated due to their 
pion-producing interactions with the CMB, i.e. the so-called Greisen-Zatsepin- 
Kuzmin cut-off [214, 215]). Pions trigger electromagnetic cascades, effectively 
transferring energy from the CRs to the GeV-TeV energy range. This can poten¬ 
tially contribute to the DGRB [216, 42, 43, 217]. Ref. [42] estimates this emission, 
showing that it is subject to large uncertainties. Indeed, the signal strongly de¬ 
pends on the evolution of the UHECRs with redshift, on their composition and 
on the intensity of the magnetic fields encountered during the propagation of the 
CRs. The uncertainty on the intensity of the emission spans more than two or¬ 
ders of magnitude below 10 GeV. The most optimistic prediction indicates that 
UHECRs can indeed represent a significant contribution to the EGRET DGRB in 
Refs. [6, 7]. However, above 10 GeV, the intensity of the signal reduce signihcantly 
so that it would not be able to explain the bulk of the sub-TeV DGRB detected 
by the Fermi LAT. 

• Type la supernovae: Type la supernovae are generated in the thermonuclear 
explosions of white dwarfs near the Chandrasekhar mass [218, 219]. Gamma-ray 
emission results from the decay of the material (mainly ®®Ni) produced during 
the detonation. Refs. [220, 221, 40] compute the cumulative emission associated 
to this class of supernovae as a function of their event rate. The latter is either 
measured directly or related to the cosmic SER, p*(^), through a model of the 
delay time, i.e. the time required for a Type la supernova progenitor to become 
a supernova. The cumulative gamma-ray emission can contribute significantly to 
the DGRB only around the MeV scale (see also Ref. [222]). 

On the other hand. Ref. [41] computes the emission associated with the decay 
of neutral pions produced in the interactions of CR protons with the interstellar 
medium of the galaxy hosting the supernova. This is a scenario very similar to the 
one described in Sec. 2.2.3 for SEGs. However, in Ref. [41], the authors consider 
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CRs accelerated in shocks induced by the explosions of Type la supernovae and not 
of the core-collapse supernovae, as done in Sec. 2.2.3. The predictions of Ref. [41] 
are affected by a significant uncertainty and they span a range between less than 
10% and 100% of the DGRB reported by the Fermi LAT in Ref. [8] between 1-10 
GeV. 

• Gamma-Ray Bursts: gamma-ray bursts are very short and intense episodes 
of beamed gamma-ray emission with a bimodal SED [223, 224, 225]. As in the 
case of AGNs, the low- and high-frequency peak are associated, respectively, with 
synchrotron radiation and IC. The widely-accepted fireball internal-external shock 
model [226] allows to describe the phenomenology of the material inside the bursts 
and to predict the SED of the bursts. Refs. [227, 228] estimate the total contribu¬ 
tion of unresolved gamma-ray bursts to the DGRB by adopting the LF given in 
Ref. [229]. Results suggest that this emission can explain only a small fraction of 
the EGRET DGRB in Ref. [6], becoming negligible at energies above ~ 40 GeV. 
Similar results have been obtained in Ref. [230] which estimates the contribution 
of unresolved gamma-ray bursts to be as large as 0.1% of the EGRET DGRB in 
Ref. [6] at the GeV scale. 

• Small Solar-system bodies: our knowledge of comets populating the Oort 
Cloud is quite limited and it comes almost entirely from numerical simulations. 
We believe that more than 10^^ objects, with sizes ranging from 1 to 50 km, may 
populate that region. These small Solar-system bodies emit gamma rays from the 
hadronic interactions of CRs impinging on them. Their abundance is the main 
unknown, with column densities of Solar-system bodies spanning over three orders 
of magnitude. As a consequence, a similar level of uncertainty affects the predic¬ 
tions for their cumulative gamma-ray emission. Their contribution to the DGRB 
may go from overpredicting the DGRB measured in Ref. [7] to being responsible 
for only few percents of it. See Ref. [231] and references therein. 

• Radio-quiet AGNs: AGNs at sub-Eddington luminosities are characterized by 
radiatively inefficient accretions. Since the jet is not beamed enough to trigger 
non-thermal emission, the source emits mainly between IR and X-ray frequencies. 
Gamma rays can still be produced: one possibility is to have proton-proton in¬ 
teractions producing neutral pions in the hot gas surrounding the supermassive 
Black Hole. The intensity of the signal depends on the spin of the Black Hole, 
since more rapidly rotating objects correspond to larger X-to-gamma flux ratio 
[232, 233, 234]. Another possibility is a population of non-thermal electrons that 
can be accelerated in the hot corona around the AGN [235, 16]. The similarities 
between solar coronae and accretion disks [236] suggest that magnetic reconnec¬ 
tion could be responsible for this acceleration [237]. The scenario is considered 
in Refs. [235, 16] where, using a luminosity-dependent-density evolution, the au¬ 
thors determine that this contribution can explain the whole DGRB measured by 
EGRET in Refs. [6, 7] , but only below 1 GeV. 

• Imperfect knowfledge of the Galactic foreground: the authors of Ref. [238] 
reconsidered the measurement of the DGRB by EGRET and noted that some 
modifications on the modeling of Galactic GRs could decrease significantly the in¬ 
tensity of the DGRB. They also raised some doubts on the approach of template 
fitting. In particular, the effect of an extended halo of electrons around the MW 
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Figure 9: The energy spectrum of the DGRB (black points) as recently measured by the Fermi LAT 
[9]. Gray boxes around each data point denote the uncertainty associated with the Galactic diffuse 
emission. The solid color lines indicate the expected gamma-ray emission from unresolved sources, for 
4 different well-established astrophysical populations: blazars (in orange), MAGNs (in green), SFGs (in 
blue) and MSPs (in red). Color bands represent the corresponding uncertainties on the emission of each 
population. Estimates are taken from Ref. [25] (blazars). Ref. [29] (MAGNs), Ref. [161] (SFGs) and 
Ref. [38] (MSPs). 


(with a consequent IC gamma-ray emission extending to high latitudes) is con¬ 
sidered. Furthermore, Ref. [239] investigates the possibility of a gas cloud with a 
mass of few IO^^Mq, extending to hundreds of kpc from the center of the MW. 
This halo would be theoretically well motivated, as it would alleviate the problem 
of the missing baryons in spiral galaxies. A similar object around spiral galaxy 
NGC 1961 would also explain the diffuse X-ray detected in Ref. [240]. Hints of 
such large halo could be already present in hydrodynamical A^-body simulations of 
our Galaxy [241, 242, 239]. The gamma-ray emission associated with pion decay 
in this hypothetical gas halo would be able to explain between 3% and 10% of the 
Fermi LAT DGRB in Ref. [8], depending on the exact size of the halo. 


Other possibilities not considered in the list above include emission from massive 
black holes at z ~ 100 [243], from the evaporation of primordial black holes [244, 245], 
from the annihilations at the boundaries of cosmic matter and anti-matter domains [246] 
and from the decays of Higgs or gauge bosons produced from cosmic topological defects 
[247], 

We conclude this section by discussing Fig. 9. The image gathers the most recent 
predictions for the “guaranteed” components to the DGRB, i.e. the emission associated 
with unresolved blazars, MAGNs, SFGs and MSPs (see sections from 2.2.1 to 2.2.4). 
They are taken from the results of Refs. [25, 29, 161, 38], respectively and they are 
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depicted in Fig. 9 by orange, green, blue and red lines, respectively.^^ Each contribution 
is embedded in a band that denotes the level of uncertainty affecting the prediction. The 
largest is the one associated with MAGNs (light green band) spanning almost one order of 
magnitude. Black data points represent the new Fermi LAT measurement of the DGRB 
in Ref. [9] (see Sec. 2.1). The gray boxes around the data points indicate the systematic 
error associated with the modeling of the Galactic foreground. From the figure, it is 
clear that MSPs are subdominant and that the remaining 3 astrophysical components can 
potentially explain the whole DGRB, leaving very little room for additional contributions 
(see also Refs. [61, 248, 217]). Similar results have been recently obtained by Ref. [65]. 
This reference also shows that the goodness of the fit to the Fermi LAT DGRB energy 
spectrum in terms of astrophysical sources depends significantly on the model adopted 
for the diffuse Galactic foreground and on the slope of the energy spectrum of unresolved 
SFGs. In particular, a description of SFGs with a softer energy spectrum (similar to 
that of the Galactic foreground) can provide a better fit to the DGRB intensity. 

2.3. The Dark Matter component of the Diffuse Gamma-Ray Background 

The DGRB can also be used to investigate more exotic scenarios than those presented 
in the above subsections. In particular, it has already been shown that the DGRB is a 
powerful tool to investigate the nature of DM. 

Discussing the very wide range of viable DM candidates is beyond the scope of this 
review (see, e.g., Ref. [249]). In the following, we only consider a family of candidates 
called Weakly Interacting Massive Particles (WIMPs), loosely characterized by a mass 
of the order of the GeV-TeV and by weak-scale interactions. This is a very well studied 
scenario since many extensions of the Standard Model of Particle Physics predict the 
existence of WIMPs [250, 251, 44, 252, 253]. It is also quite natural for WIMPs to 
reproduce the DM relic density observed, e.g., by Planck [254]. Yet, currently there is 
no observational confirmation of the existence of WIMPs. 

WIMP DM can either annihilate or decay into Standard Model particles, including 
gamma rays. This is a general prediction of WIMP candidates and it represents an 
additional reason to focus only on WIMPs for this review. The specific mechanisms 
of gamma-ray emission (see, e.g.. Ref. [44] for a review) depend on the DM candidate 
considered and include i) direct production of monochromatic gamma rays, ii) decay of 
neutral pions, produced by the hadronization of the primary annihilation/decay prod¬ 
ucts, Hi) final state radiation and iv) secondary emission by IG or bremsstrahlung of 
primarily produced leptons. Since no DM source has been unambiguously detected up 
to now, the entire DM-induced gamma-ray emission may be unresolved and, thus, it 
contributes to the DGRB. In Sec. 2.3.1 we discuss the potential DM contribution to 
the DGRB in the case of self-annihilating DM particles, while Sec. 2.3.2 is devoted to 
decaying DM. Note that some DM candidates can experience both annihilations and 
decays [255]. 

DM-induced gamma rays can be produced in the DM halo of the MW or in ex- 
tragalactic DM structures and substructures. We refer to the two possibilities as the 
“Galactic” and “cosmological” DM components, respectively. The latter is isotropic by 


'^^Ref. [25] only provides the total emission from resolved and unresolved blazars. Since we are inter¬ 
ested in the unresolved component, the orange line in Fig. 9 is obtained by subtracting the emission of 
resolved sources from Ref. [9] from the total sigual from blazars. The width of the light orange band is, 
then, computed summing the estimated errors of the two components in quadrature. Also, the abrupt 
end of the SFG component at ~ 100 GeV in not physical and simply comes from the lack of predictions 
above this enegy in Ref. [161] only up to that energy. 
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construction, while the former is expected to exhibit some anisotropy, due to the partic¬ 
ular location of the Earth in the DM halo of the MW. We remind that, as described in 
Sec. 2.1, the intensity of the DGRB is obtained by means of an isotropic template [9]. 
However, the Galactic DM signal can exhibit a significant anisotropy and, in that case, 
it cannot be considered part of the DGRB.^® In this section, we focus mainly on the 
contribution of the cosmological DM signal to the DGRB, discussing a possible Galactic 
DM component only towards the end of the section. 


2.3.1. The case of annihilating Dark Matter 

In the ACDM cosmological framework [254] , initial matter fluctuations in the Early 
Universe are the seeds of the structures that populate today’s Universe. These fluctu¬ 
ations grow by accreting new matter and form the first protostructures, which, then, 
collapse and eventually virialize into DM halos. ACDM predicts that, in later epochs, 
larger halos gradually assemble by accretion and merging of smaller halos. Under this 
scenario of structure formation, a cosmological gamma-ray signal is expected from the 
annihilations of DM particles taking place in all DM halos at all cosmic epochs. 

The cosmological gamma-ray flux d^/dEodkl (i.e. the number of photons per unit 
energy, time, area and solid angle) produced by DM annihilations at energy Eq over all 
redshifts z is given by [48, 49, 51]: 


dE^dO. 
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(15) 

where is the mass of the DM particle, {av) its annihilation cross section^^ and the 
sum runs over all the possible annihilation channels, each of them corresponding to a 
specific branching ratio Bi and a differential photon yield dN^/dE. is the current 
DM density ratio, pc the critical density of the Universe, and H(z) and c are the Hubble 
parameter and the speed of light, respectively. The function tebl(R^) z) accounts for the 
absorption of gamma-ray photons due to interactions with the EBL. Finally, the quantity 
C(z) is the so-called flux multiplier and it indicates the variance of the fluctuations in 
the field of squared DM density. It is, therefore, a measure of the statistical clustering 
of DM in the Universe: 

9 {P%{^)) 


where is the comoving DM density and the parentheses (•) denote angular integration 
over all the possible pointings in the sky. 

Two approaches have been proposed to calculate the flux multiplier. The first is 
based on the Halo Model [200, 201] and it relies on the knowledge of the abundance 
and of the internal properties of DM halos. The former is described by the Halo Mass 
Function (HMF), dn/dM, i.e. the comoving number of halos per unit mass, while the 
latter is encoded in the DM halo density profile. More specifically, in the Halo Model 
framework, (^{z) is proportional to the integral of the HMF times the integral of the DM 


Small-scale anisotropies in the DGRB, on the other hand, are discussed in detail in Sec. 3. 

^^The annihilation cross section in Eq. (15) is assumed to be dominated by s-wave interactions. In 
the case of a dependence on the relative velocity of the annihilating DM particles, Eq. (15) has to be 
modified accordingly and the signal will, thus, depend on the velocity distribution of DM [256, 257]. 
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density squared inside the halo [48, 49, 50, 64]: 


dM^M^{F(M,z)), (17) 

Pc dM 3 

where Madn is the minimum halo mass. Typical values for range approximately 

between 10“^^ and 10 “^Mq for supersymmetric DM candidates.^® Its exact value de¬ 
pends on the position of a cut-off in the power spectrum of matter fluctuations, above 
which the formation of DM structures is suppressed. This cut-off arises from the com¬ 
bined effect of kinetic decoupling and baryonic acoustic oscillations [258, 259, 260, 261], 
and its precise location ultimately depends on the Particle Physics nature of the DM 
candidate [262, 261, 263]. In addition, the mapping between the cut-off scale and Mmin 
is not well dehned, depending on the assumed relation between mass and size of the 
small-mass halos [264]. These uncertainties are responsible for the huge variability of 
Mmin quoted above. 

The factor A{z) in Eq. (17) is the so-called halo overdensity and it depends on the 
cosmology assumed and on the details of the gravitational collapse of the halo. The 
radius Ra at which the mean enclosed density of a DM halo is A( 2 ;) times pc is called 
the virial radius, which can be taken as a measurement of the size of the halo. The DM 
distribution inside the halo is codified in the function {F(M, z)): 


F{M,z) = cl{M,z) 


dxx^K?{x) 
dx k{x)]‘^ 


(18) 


where k{x) is the DM density prohle and x = r/r^. is a scale radius, whose precise 
definition depends on the assumed k{x). The quantity ca is the so-called concentration 
[265, 266, 267] and it is defined as Ra/^s- Note that {F{M, z)) in Eq. (17) is the function 
F{M,z) from Eq. (18) averaged over a log-normal probability distribution assumed for 
CA- Such a distribution accounts for halo-to-halo scatter on the value of ca [265, 268, 
269], which is a natural consequence of the stochastic process of structure formation in 
ACDM cosmology. 

Information on the abundance and internal properties of DM halos is mainly ex¬ 
tracted from iV-body cosmological simulations (see, e.g., Ref. [270] and references therein) 
However, simulations do not resolve the whole halo hierarchy down to Mmin- The current 
resolution limit for MW-like DM halos is approximately at IO^Mq at 2 ; = 0 [271], i.e. 
several orders of magnitude above the expected value for Mmin- Thus, extrapolations 
of the relevant quantities (e.g., the HME and ca(M, z)) are needed. Since F(M,z) in 
Eq. (18) depends on the third power of the halo concentration, the way ca(M, z) is 
extended below the mass resolution of simulations is crucial in the estimation of the cos¬ 
mological DM signal. Above the mass resolution limit, ca(M) can be well described by a 
power law [265, 272, 273, 266, 267]. Several works assume the same behavior to be valid 
down to Mmin [274, 56, 275, 276], thus assigning very large concentrations to the smallest 
halos. This translates into a very high gamma-ray flux expected from DM annihilations. 
Nevertheless, power-law extrapolations for ca(M) are not physically motivated and they 
are now clearly ruled out both by recent high-resolution simulations of the smallest DM 
halos [277, 278, 279] and by theoretical predictions deeply rooted in the ACDM cosmo¬ 
logical framework [267, 280, 281]. Indeed, simulation-based and theoretical estimates of 


^®The value Mmin = 10 ^Mq has become a standard benchmark value in the field. 
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ca{M, z) have been shown to agree now remarkably well over the full halo mass range, 
i.e. from Earth-like “microhalos” up to the scale of the heaviest galaxy clusters [280]. 
Compared to power-law extrapolations, these estimates exhibit a flattening of c^{M) at 
low masses. This leads to substantially less concentrated low-mass halos and, thus, to a 
considerably smaller cosmological DM signal. Overall, predictions for the cosmological 
DM-induced emission can vary by few orders of magnitude depending on the adopted 
model for ca{M,z) in Eq. (18), see, e.g.. Refs. [56, 282]. 

Not unexpectedly, the variability of the predicted DM signal also depend on the 
particular choice of M^\^. Eor example, the cosmological signal increases by up to a 
factor of ~ 6 when Mmin goes from 10 “^Mq//i to 10 “^^Mq//i [62]. This refers to the 
case of a power-law extrapolation of ca(M) and, by construction, flux multipliers that 
rely on power-law extrapolations are particularly sensitive to the choice of Mmin- When 
adopting a ca(M) that flattens at low halo masses, the flux multiplier changes by a 
factor of ~ 3 over the same range of Mmin [62, 64]. 

A^-body cosmological simulations have also been employed to understand the HME 
and its redshift evolution [283, 284]. Mock halo catalogs have been used to test the 
predictions of the Press-Schechter formalism, according to which the HME can be written 
as follows [200, 285, 286, 287]: 


dn 

dM 


{M,z) 


f{a{M,z)) 


Pc0.y^{z) dlna ^(M, z) 
M dM 


(19) 


where cr{M, z) is the variance of the fluctuations of the DM density field (smoothed on a 
scale of mass M) and the exact expression for f{x) depends on the mechanism assumed 
to describe the halo gravitational collapse. An accurate htting formula for the HME, 
inspired by Ref. [288], can be found in Ref. [289] for the cosmological model favored 
by the hrst data release of the Wilkinson Microwave Anisotropy Probe (WMAP). More 
recently. Ref. [267] have also adopted the functional form proposed in Ref. [289] and 
derived the parameters of the HME compatible with the cosmological model preferred 
by Planck.^® Overall, the HME at z = 0 can be qualitatively approximated by a power 
law with a slope between -1.9 and -2.0 and a sharp cut-off for halos more massive than 
~ IO^^Mq. The uncertainty in the calculation of the cosmological DM signal induced 
by different possible parametrizations of the HME is only marginal when compared to 
other sources of uncertainties. Eor instance, assuming the HME of Ref. [288] instead of 
the one in Ref. [289] changes the total gamma-ray flux by a factor of about 20% (see 
also Ref. [282]). 

High-resolution A^-body simulations have also helped to establish the density profiles 
of DM halos in great detail. Ref. [291] determined that DM halos exhibit a density 
profile that is universal, i.e. independent of the halo mass. A Navarro-Erenk-White 
(NEW) profile, i.e. proportional to k oc 1/r in the inner region and to a steeper r~^ at 
large radii, provides a good ht. More recently, the so-called Einasto profile was found to 
agree better with the results of A^-body simulations, especially at intermediate halo radii 
[292, 293, 294]. Even if many other parametrizations have been proposed over the last 
years [295, 296, 297, 298, 299, 300], current A'-body simulations agree on a slope < —1 
for the DM density in the inner region. However, these “cuspy” profiles are derived from 
simulations that only contain DM, without including baryons. Indeed, high-resolution 


^®Note that the agreement of different fitting formulas with the simulations may depend on the algo¬ 
rithm used by the simulators to extract the mass of a halo from the raw particle data, with the caveat 
that different algorithms may lead to different results [290]. 
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observations of the rotation curves of DM-dominated dwarfs and low-surface-brightness 
galaxies favor DM density profiles with a flat central core [301, 302, 303, 304, 305, 
306, 307, 308, 309]. Phenomenological cored profiles, such as the Burkert one [302], were 
proposed to accommodate such results. More recently, hydrodynamical simulations have 
been performed, which realistically include the complexity of baryonic physics. They 
begin to reproduce the observed properties of galaxies successfully, e.g. in Refs. [310, 
311]. Yet, the exact interplay between baryons and DM at all radial scales and for all 
halo masses is not fully understood, and the impact of complex baryonic phenomena 
(such as supernova feedback, stellar winds and baryonic adiabatic compression) on the 
DM density profile is still uncertain, particularly in the inner region [312, 313, 314, 315, 
316, 317, 318, 319, 320]. Ref. [321] finds a difference of almost one order of magnitude 
in the cosmological flux when a Burkert profile is assumed for all DM halos instead of a 
NFW one. 

In addition to DM halos, a natural prediction of ACDM is the existence of a large 
number of subhalos, i.e. halos gravitationally bound to a larger host halo and located 
within its virial radius. Since the annihilation luminosity of a halo is proportional to 
its DM density squared, the presence of small clumps with high DM densities has the 
effect of boosting the overall gamma-ray luminosity of the host halo [322, 323, 324, 274, 
325, 326, 327, 56, 328, 329, 275, 330, 276, 279, 331, 280]. The additional contribution 
from substructures can be accounted for in the computation of the flux multiplier by 
adding an extra term in Eq. (17). In particular, the factor {F{M, z)) has to be replaced 
by {F{M^z)){l + B{M,z)), where B{M,z) is called the “boost factor”. 

The subhalo population is characterized by a subhalo mass function: 

extending down to Mmin- The slope of the subhalo mass function has been measured in 
high-resolution A^-body simulations, ranging between -1.9 and -2 [332, 271]. These values 
are also in line with theoretical expectations from the Press-Schechter theory of structure 
formation [333, 334, 335]. It has been noted, though, that several effects may prevent the 
subhalo mass function to reach the lowest subhalo masses, since processes such as tidal 
disruption, accretion or merging may be particularly efficient and deplete the low-mass 
tail of the mass function. It is difficult to estimate the survival probability of these small 
subhalos, while it remains computationally very expensive to simulate and keep track 
of such processes with the resolution needed. Although the properties of low-mass DM 
subhalos are expected to follow those of the more massive counterparts, the abundance 
and distribution of DM substructures below the resolution of current simulations remain 
uncertain. The internal properties of subhalos are also subject to debate. We still lack a 
comprehensive understanding of the subhalo concentration, even if subhalos have been 
shown to exhibit larger C/\{M) than halos of the same mass [336, 271, 337, 338]. Similarly 
to the case of main halos, the assumptions made on ca(M) and on the abundance of low- 
mass subhalos have a large impact on the amplitude of the subhalo boost factor. Overall, 
under the most extreme scenarios (corresponding to blind power-law extrapolations of 
ca(M)), the presence of subhalos can increase the total cosmological annihilation signal 


^°This definition of the boost factor is particularly convenient for the computation of the flux multiplier 
and of the cosmological DM signal. We warn the reader, though, that alternative definitions of B{M, z) 
can be found in the literature. 
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by more than one order of magnitude [48, 62, 217, 339]. 

As an alternative to the Halo Model described above, the Power Spectrum approach 
has been recently introduced to compute the flux multiplier [340, 264]. In this new 
framework, ((z) can be calculated by means of the non-linear matter power spectrum 
PjVL, he. the Fourier transform of the two-point correlation function of the matter 
density field, as follows: 


C(^) = (5^(^)) 



dk k^P^hik, z) 
k 27r^ ’ 


( 21 ) 


where femax('2^) is a maximal scale that can be related to in the Halo Model formal¬ 

ism.^^ The main benefit of the Power Spectrum approach relies on the fact that only one 
single quantity is needed for the calculation of C(-2), i-e. Pnl in Eq- (21). Pnl can be 
measured directly in A^-body cosmological simulations using only a matter density map 
and it does not rely on the concept of DM halos. It is, thus, also independent on the 
complex issue of halo finding and on the uncertainties associated with it [341].^^ This 
is, indeed, what the authors of Ref. [264] did, guided by the results of the Millennium-I 
and Millennium-H simulations [284, 343]. Yet, as in the Halo Model approach, current 
A^-body simulations do not reach the highest k values needed in Eq. (21) and, thus, ex¬ 
trapolations beyond the simulation resolution are again required.Note, however, that 
only Pnl is extrapolated in this case, as opposed to the Halo Model approach that relies 
on the knowledge of several quantities in the low-mass regime. This is indeed the reason 
why the authors in Ref. [64] adopted the Power Spectrum approach to obtain a realistic 
estimate of the uncertainty on C(^)) even if their fiducial flux multiplier is computed 
within the Halo Model framework. Furthermore, within the PS approach, it is possible 
to motivate in a realistic way both the choice of fcmax and the way the extrapolation 
is performed. Note also that, by construction, the Power Spectrum approach naturally 
includes the contribution of substructures down to length scales ~ vr/fcmax, while in the 
Halo Model the description of snbstructures requires additional knowledge, often leading 
to further debatable extrapolations and uncertainties.^^ 

A comparison between the Power Spectrum and Halo Model approaches is performed 
in Refs. [264, 344, 64]. Results from Refs. [264, 64] are reproduced in Fig. 10. In the left 
panel, the red bands indicate the predictions for Q{z) (multiplied by factors depending 
on redshift) obtained by means of the Power Spectrum method for different choices of 
kma.x- For each value of /cmax adopted, the corresponding band accounts for different 
ways of extrapolating Pnl beyond the resolution of the Millennium simulations. The 
predictions for the Halo Model (gray band) refer to a Mmin = 10“®Mq//i and have been 
obtained following Refs. [56, 58]. In the right panel of Fig. 10, the predictions of the 


^^The exact relation between and Mmin is not trivial, as discussed in detail in Refs. [264, 64]. 

^^In Refs. [342, 331], the authors introduce the concept of particle phase space average density, an 
estimate of the coarse-grained phase-space density of DM structures. Interestingly, they show in Ref. [331] 
how this observable can be used to determine the DM-induced emission of MW-like halos. The particle 
phase space average density can be directly computed directly from A^-body simulations’ raw data and, 
therefore, the formalism shares some similarities with the Power Spectrum approach, introduced here to 
estimate the cosmological DM signal. 

^®At redshift zero, for instance, the maximum k values resolved in current large-scale structure simu¬ 
lations are typically of the order of few hundreds, while fcmax may take values up to 10®. 

^^We remind the reader that Eq. (21) only determines the cosmological DM signal. Any gamma-ray 
emission associated with the halo of the MW and its subhalos is not accounted for by the Power Spectrum 
model and, thus, needs to be added by hand (see further discussion at the end of this section). 
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Figure 10: Comparison of the Halo Model and Power Spectrum approaches for the calculation of the 
flux multiplier. The Halo Model makes use of Eq. (18) while the Power Spectrum formalisms adopts 
Eq. (21). Left: Power Spectrum predictions are shown in red. The different bands correspond to different 
prescriptions to derive the cut-off scale fcmax from the size of the DM halo with mass The width 

of the bands accounts for different extrapolation schemes up to fcmax. The Halo Model predictions 
(gray band) are obtained following Refs. [56] and [58]. Figure taken from Ref. [264]. Right: Predictions 
according to the Power Spectrum approach (gray) and the Halo Model (red). The Halo Model formalism 
follows what done in Ref. [64] for two possible subhalo mass functions. The gray band corresponds to 
a particular choice of fcmax in the Power Spectrum approach, for different ways (labeled as “PS (max)” 
and “PS (min)” in the figure) to extend Pnl beyond the resolution of the simulations. For both panels, 
a value of Mmin = 10~^Mq/H is adopted in the Halo Model formalism. Figure taken from Ref. [64]. 


Power Spectrum approach (for a specific value of fcmax) are given by the gray band. The 
width of the band corresponds to two different prescriptions to perform the extrapolation 
(labeled “PS (min)” and “PS (max)” in the figure). On the other hand, the red band 
reproduces the results of the Halo Model in Ref. [64] for Mmin = W~^MQ/h and two 
values of the slope for the subhalo mass function. Remarkably, the two methods agree 
well in their predictions for C('2^)) within uncertainties. This is so despite the caveats 
mentioned above and the intrinsic differences of the two formalisms. 

From Sec. 2.2 we know that gamma-ray emission of astrophysical origin is able to 
explain a significant fraction, if not all, of the DGRB, therefore leaving very little room 
for a potential DM contribution (see also Fig. 9). Moreover, the gamma-ray energy 
spectrum expected from DM annihilations is generally slightly curved, with a cut-off 
at the mass of the DM particle and the possibility of spectral features like bumps or 
lines. On the other hand, the DGRB detected by the Fermi LAT exhibits a power- 
law spectrum at low energies and a cut-off compatible with EBL attenuation at higher 
energies (see Sec. 2.1). This suggests, again, that DM annihilations cannot provide a 
dominant contribution to the DGRB. Thus, the DGRB can be used to set limits on the 
intensity of the DM-induced emission. These are usually translated into upper limits 
in the (rriy-, {crv)) plane, identifying the region that, given a model for the abundance 
and properties of DM halos and subhalos, is excluded as it overproduces the measured 
DGRB. 

Most of the works following this idea employ the Halo Model approach to predict the 
cosmological DM annihilation signal. Given the large number of parameters involved in 
the Halo Model framework, as well as the large uncertainties associated with some of 
them, it is very hard to perform a detailed, one-to-one comparison among the different 
DM limits available in the literature. In particular, predictions obtained by different 


35 




















groups differ mainly due to different assumptions on ca{M,z) and Mmin - In addition, 
some works also consider the Galactic DM signal. Differences can be further amplihed 
by the various statistical prescriptions employed to compute the upper limits. Yet, we 
believe that a comparison among the limits in the literature is useful, as it showcases 
the potential of searching for DM in the DGRB compared to other indirect DM probes. 
Fig. 11 summarizes some of the upper limits available.^® They all refer to annihilations 
entirely into b quarks. Note that the limits obtained by assuming a power-law ca(M) 
below the mass resolution of Al-body simulations are among the most constraining in 
Fig. 11 (see, e.g., the gray dashed line from Ref. [217]). However, as argued above, power- 
law ca{M) are not well motivated, putting the corresponding limits into question. 

Two approaches are possible when deriving DM limits from the DGRB measurement. 
In the first one it is required that the DM signal does not overshoot the measured DGRB 
at a particular confidence level, typically 2 or 3a. See, e.g.. Refs. [365, 59, 366, 60, 62, 335, 
339]. This leads to conservative and robust upper limits on the DM annihilation cross 
section, shown by the green lines in Fig. 11. Alternatively, one or more astrophysical 
contributions to the DGRB can be modeled and included in the analysis. The DM 
component is, then, required not to overshoot the fraction of the DGRB not already 
accounted for by astrophysics. These limits are more constraining than the previous ones 
[58, 367, 61, 109, 248, 217]. Nevertheless, since the exact contribution from astrophysical 
sources is not known, it must be noted that they are subject to larger uncertainties. 

The most recent constraints on annihilating DM comes from Ref. [64] and are shown 
by black lines in Fig. 11. The “conservative limit” (solid black line) is obtained without 
including any modeling of astrophysical contributors to the DGRB. These limits exclude 
annihilation cross sections a factor of ~ 3 lower than the thermal value of 3x 10“^®cm^s“^ 
for a mass of 10 GeV, while, for > 1 TeV, the upper limit is around 5 x 10“^^cm^s“^. 
Ref. [64] also shows that an improvement of one order of magnitude (a factor ~ 2) 
for DM masses of the GeV scale (around 30 TeV) is possible when i) the cumulative 
astrophysical contribution to the DGRB is modeled as a power-law in energy with an 
exponential cut-off and ii) slope and position of the cut-off are fixed to the best-fit values 
to the DGRB data in Ref. [9]. Note that this approach (referred to as “sensitivity reach” 
and denoted by the dashed black line in Fig. 11) does not account for any uncertainty 
in the description of the astrophysical emission. A more realistic scenario is the one 
of Ref. [25], where the authors estimate the unresolved astrophysical emitters using 
the most up-to-date information from resolved sources or from other frequencies (see 
Sec. 2.2). A renormalization factor A is included in front of the total astrophysical 
contribution to account for possible fluctuations in its intensity. The limit of {av) (solid 
blue line in Fig. 11) is then obtained by profiling over A. This results in an improvement 
of a factor of ~ 3 for niy. = 10 GeV, with respect to the conservative limits of Ref. [64]. 
A negligible improvement is expected at masses larger than 10 TeV, where the limits 
are determined by the Galactic DM signal. 

Similar results have recently been obtained in Ref. [65], which also employs a model 
for the astrophysical components to the DGRB. The different statistical analysis per¬ 
formed in Ref. [65] suggests that a description of the DGRB of Ref. [9] in terms of 
unresolved astrophysical sources and gamma rays from DM annihilations in the MW 


only consider limits derived from the two measurements of the DGRB performed by the Fermi 
LAT [8, 9]. Older works based on earlier DGRB measurements can be found, e.g., in Refs. [345, 346, 
347, 48, 49, 348, 349, 350, 50, 351, 51, 52, 352, 353, 354, 355, 356, 274, 357, 358, 55, 359, 321, 360, 361, 
362, 363, 56, 364]. 
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Figure 11: Upper limits obtained by considering the DGRB energy spectrum measured in Refs. [8, 9]. 
Annihilations into b quarks are assumed. The regions above the colored lines are excluded because the 
cumulative DM-induced emission would overproduce the DGRB. Different lines correspond to different 
assumptions for the properties of DM halos (especially for low halo masses) and different methods to 
compute the upper limits. The solid green line is taken from Fig. 2 of Ref. [365] while the dashed green 
line is from Fig. 8 of Ref. [335] (lower bound of the band relative to am = 2 for the emission from the 
Galactic Poles). The dashed blue line is taken from Fig. 5 of Ref. [58] (conservative limits for model 
MSII-Subl) and the solid gray one from Fig. 2 of Ref. [367]. The solid and dashed red lines are taken 
from Fig. 3 of Ref. [61] (limits labeled “Fermi EGB”) and from Fig. 5 of Ref. [248] (panel labeled “best-fit 
background”), respectively. The dashed gray line is from Fig. 15 of Ref. [217] (default substructures’ 
model). The black lines correspond to the predictions obtained in Ref. [64] by means of the Halo Model 
(reference scenario). The solid blue line is taken from Fig. 4 of Ref. [25], while the dashed purple one is 
from Fig. 4 of Ref. [65]. The blue region indicates the portion of the {m^, {crv}) plane already excluded by 
the observation of the Segue 1 dwarf Spheroidal galaxy performed by the MAGIC telescopes (see Fig. 6 
of Ref. [368]). The dark gray region is excluded by the analysis performed by the H.E.S.S. telescopes in 
Ref. [369] from the so-called “Galactic Center halo” (see their Fig. 4 for an Einasto DM density profile). 
Finally, the light gray region indicates the DM candidates not compatible with the combined analysis of 
15 dwarf Spheroidal galaxies by the Fermi LAT [370]. A comparison between the Fermi LAT DGRB and 
the DM-induced signal can also be found in Refs. [366, 59, 60, 109, 62, 339]. The dash-dotted horizontal 
line marks the value of the thermally averaged annihilation cross section. 
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halo provides a better fit, compared to a purely astrophysical interpretation with no 
DM. Such a hint of a DM signal in the DGRB suggests a DM particle with a mass of 
~ 10 — 20 GeV (for annihilation into b quarks), depending on the model adopted for the 
diffuse Galactic foreground.^® 

Indeed, some of the limits shown in Fig. 11 have been obtained assuming that DM 
annihilations in the MW halo also contribute to the DGRB [58, 365, 367, 59, 366, 61, 109, 
62, 335, 339, 248, 217]. Note that, as mentioned in Sec. 2.1, the DGRB is obtained from 
the normalization of the isotropic template in the multi-component ht to the gamma- 
ray data at high Galactic latitudes. Thus, when a Galactic DM signal is included, it is 
implicitly assumed that the DM signal is sufficiently isotropic in this particular region 
of the sky. 

Two distinct components contribute to the Galactic DM signal. The first accounts 
for the smooth DM distribution of the host DM halo of the MW, while the second one 
comes from the population of Galactic subhalos. The former depends on the DM density 
profile assumed for the MW main halo. This is uncertain in the innermost parts of the 
Galaxy, but for |6| > 20°, i.e. more than 3 kpc from the Galactic Center, iV-body 
simulations roughly agree. The signal from this smooth component peaks towards the 
Galactic Center and, between b = 20° and 90°, typical variations are of a factor of 
~ 16 [64]. Thus, even outside the Galactic plane the signal is largely anisotropic and, 
therefore, this component cannot be described by an isotropic template as the DGRB. 
Indeed, the presence of an emission with such a well-defined morphology may impact 
the procedure used in Ref. [9] to measure the DGRB energy spectrum. This was tested 
in Ref. [64], where the authors re-derived the DGRB including an additional template 
for the smooth Galactic DM signal. They found that this signal can be degenerate 
with other diffuse Galactic emissions, especially the one from IC. They also checked 
the impact that the new template would have on the upper limits on (<tu). The result 
suggests that, at least for DM candidates not excluded by the conservative limits in 
Ref. [64], this additional Galactic DM template has only a moderate effect. 

The second contribution to the Galactic DM signal comes from the subhalos of the 
MW: the brightest or closest of them may potentially give rise to bright spots in the 
gamma-ray sky. However, none of the unassociated sources in the 2FGL catalog has 
been robustly interpreted as a DM subhalo [371, 372]. The overall subhalo population 
is expected to give rise to a diffuse smooth emission [274, 53, 55, 359]. Its morphology 
depends on the abundance and distribution of subhalos in the Galaxy. A general predic¬ 
tion is that the cumulative emission of subhalos is more extended (thus more isotropic) 
than that of the main halo. More specifically, factors between 0.1 and 2 are quoted in 
Ref. [64] for the variation of the signal of Galactic substructures between b = 20° and 
90°. These small variations motivate the authors of Ref. [64] to assume this component 
as isotropic and, thus, to include it when setting their DM limits. Its impact can be 
quite signihcant since Galactic subhalos can boost the Galactic DM signal by a factor 
of 3 to 15, depending on the slope of the subhalo mass function [280]. 

We end this section by comparing the upper limits on the annihilation cross section 
derived in Ref. [64] to the results of other indirect searches for DM. In particular, in 
Fig. 11, the shaded regions indicate which portion of the parameter space has been 
already excluded by these other probes. The blue region is derived from the observation 
of the Segue 1 dwarf Spheroidal galaxy by the MAGIC telescopes [368] , while the dark 


^®A similar indication of a DM component to the DGRB was also reported in Refs. [348, 357] based 
on the EGRET measurement of the DGRB in Ref. [7] . 
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gray region indicates the portion of the (m^, {av)) space not compatible with the analysis 
of H.E.S.S. data from the so-called “Galactic Center halo” [369]. Finally, the light gray 
area is excluded by the non-detection of gamma rays from the observation of 15 dwarf 
Spheroidal galaxies with the Fermi LAT [370]. Note that the conservative upper limits 
derived in Ref. [64] (solid black line) are always inside the area already excluded, while 
the most stringent sensitivity reach (dashed black line) provides the strongest constraints 
on {av) for DM masses up to 1 TeV, above which the limit from H.E.S.S. becomes more 
stringent. 


2.3.2. The case of decaying Dark Matter 

Decaying particles can be a viable DM candidates if their decay lifetime is larger than 
the age of the Universe [373, 374, 375] . As in the previous section, we will not discuss here 
the models that can accommodate such particles or the mechanisms that guarantee long 
lifetimes. We simply consider generic WIMPs that decay emitting gamma-ray photons. 

In contrast to the case of annihilating DM, the gamma-ray signal expected from 
decaying DM particles is linearly proportional to the DM density. For instance, the 
contribution of the MW halo (at a certain energy Eq and towards the direction ip) can 
be written as follows: 


dp 

dE^dFl 


1 1 
47r 


ds PMw[r{s,ip)] E 


dN^ 


dE 


E=Eq 


( 22 ) 


where r is the DM particle lifetime, pMW the DM density profile of the MW halo and s 
is the line-of-sight variable pointing towards the direction of observation. In the case of 
prompt emission, the s-wave annihilation of 2 non-relativistic DM particles with a mass 

is equivalent to the decay of a s-wave state with mass 2m^ and integer spin [376]. 
Under those assumptions, the photon yield dN^'/dE in Eq. (22) can be determined from 
the same quantity in Eq. (15). A decaying DM candidate can also produce gamma rays 
through semileptonic channels (provided that it is characterized by a semi-integer spin) 
or through a three-body decays [255]. Mono-chromatic lines are also a typical signature 
of a vast class of decaying DM candidates [377, 378, 379, 374, 380], while the so-called 
“box-shaped” spectra have been recently introduced in Ref. [381]. Finally, secondary 
gamma rays can also be produced when primary leptons from DM decay interact via 
bremsstrahlung with the interstellar medium of the Galaxy or via IG off its radiation 
fields [382]. 

Given the linear dependence on the DM density in Eq. (22), the gamma-ray signal 
expected from the smooth MW DM halo has a different morphology, compared to the 
case of annihilating DM. For a decaying particle, the emission is now more isotropic 
and, thus, there exists a better motivation to include it among the contributors to the 
DGRB. Note also that, on the case of decaying DM, the emission of a DM halo is pro¬ 
portional to the its total mass. This means that DM substructures are not expected to 
boost significantly the predicted DM signal. Therefore, for a specific DM candidate, the 
main uncertainty affecting the emission in Eq. (22) comes from the unknown DM density 
profile of the MW halo.^”^ Remarkably, this translates into a considerably smaller uncer¬ 
tainty on the intensity of the DM-induced emission compared to the case of annihilating 
DM. 


^^When considering secondary emission, one should also add the uncertainty associated with our 
imperfect knowledge of the propagation of charged particles in the MW [382]. 
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Figure 12: Lower limits obtained by considering the DGRB energy spectrnm measured in Ref. [8]. 
Decays into are assumed. The regions below the colored lines are exclnded because the cumulative 

DM-induced emission would overproduce the DGRB. Different lines correspond to different assumptions 
for the properties of DM halos and different methods to compute the lower limits. The solid blue line is 
taken from Fig. 3 of Ref. [383] , while the solid gray one is from Fig. 8 of Ref. [384] (NFW DM density 
profile). The solid red and solid black lines are taken from Fig. 4 of Ref. [385] and from Fig. 6 of Ref. [59] 
(NFW DM density profile), respectively. The dashed gray line is from Fig. 2 of Ref. [386] (DM signal 
and power-law background) and the dashed red one from Fig. 4 of Ref. [380]. The dashed black line is 
taken from Fig. 4 of Ref. [367]. The dashed bine line comes from Fig. 3 of Ref. [387] (for model A of the 
Galactic foreground). The blue region indicates the portion of the {m^,T) plane already excluded by 
the observation of the Ursa Minor dwarf Spheroidal galaxy performed by the Fermi LAT (see Fig. 3 of 
Ref. [388], for the case of no IC emission). The dark gray region is excluded by the analysis performed 
in Ref. [380] combining Fermi LAT data from the position of 8 galaxy clusters (see their Fig. 4). Finally, 
the light gray region indicates the DM candidates not compatible with the observation of the so-called 
“Galactic Genter halo” performed by the Fermi LAT in Ref. [389] (see their Fig. 5 for a NFW DM 
density profile). 
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The cosmological signal for a decaying DM candidate is written as follows: 
d(t) _ c r ^e-"'EBL(-Eo,D cfiv 

dE^dO. Att rniDM J dE 

There is no need of including the flux multiplier since the average emission depends on 
the total amount of DM, accounted for by Ddm- This eliminates any dependence on 
T/min or on the shape of ca{M,z) and it translates into more robust predictions.^® 

Many works have used the DGRB to derive constraints on the nature of decaying DM 
candidates. Refs. [390, 391, 392, 377, 378, 393, 394, 374] consider the DGRB measured 
by EGRET in Refs. [6, 7], while Refs. [386, 385, 383, 380, 59, 367, 387] rely on the 
DGRB measurements by the Fermi LAT. Results are summarized in Eig. 12, which 
collects lower limits on r as a function of rriy^. The regions below the lines in Fig. 12 are 
excluded since the corresponding DM particle (for a specific decay channel) produces a 
gamma-ray emission which is not compatible with the DGRB data. The scatter among 
the different lines in the figure is due both to the different statistical techniques used to 
derive the limits and to the different modeling adopted for the DM-induced emission. It 
should be noted that all lines are obtained by assuming that the DM-induced emission 
is the only component of the DGRB. The only exceptions are the dashed gray line 
from Refs. [386] and the dashed blue one from Ref. [387], in which the authors also 
model the astrophysical component of the DGRB. Both the Galactic and cosmological 
signals are considered when deriving all lower limits in Fig. 12, but the ones given by 
the solid gray and dashed blue lines, which correspond to Ref. [384] and Ref. [387], 
respectively. In these two cases, limits refer to the cosmological DM component only. 
Other predictions for the contribution of decaying DM to the DGRB can be found in 
Refs. [395, 396, 62, 397]. 

The most constraining lower limit in Fig. 12 comes from Ref. [386] (dashed gray line) 
below ~ 2 TeV, and from Ref. [387] (dashed blue line) for larger DM masses. Decay 
lifetimes as large as 2 x 10^® s are excluded for DM masses below ~ 500 GeV, while the 
limit goes up to 10^® s for ~ 10 TeV. Note that the results of Ref. [387] are obtained 
from the most recent Fermi LAT measurement of the DGRB reported in Ref. [9]. The 
authors also shows how the lower limit changes depending on the model employed to 
describe the diffuse Galactic foreground. Furthermore, the lower limit at large DM 
masses is found to vary by up to a factor of a few for different ways of parametrizing 
the astrophysical component of the DGRB and its uncertainty. 

Compared to the limits derived from other DM targets, the dashed gray and dashed 
blue lines in Fig. 12 represent the most constraining information available on r. In 
particular, in Fig. 12 we show the regions excluded by three different analyses performed 
with Fermi LAT data, namely the observation of i) the Segue 1 dwarf Spheroidal galaxy 
[388] (blue region), ii) 8 galaxy clusters [380] (in dark gray) and Hi) the “Galactic Center 
halo” [389] (in light gray). 

3. The angular power spectrum of anisotropies 

In this section we move away from the study of the all-sky average of the DGRB, 
focusing on what can be learnt from its spatial fluctuations. Note that, following the 



^®Even if one rewrites Eq. (23) in terms of the Halo Model formalism, the result of integrating the 
HMF times the number of gamma rays expected from a DM halo of mass M will be quite insensitive to 
the value adopted for Mmin [62]. 
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procedure outlined in Sec. 2.1 and used in Ref. [9], the DGRB should be isotropic by con¬ 
struction. Yet, the template fitting is not sensitive to moderate small-scale anisotropies 
in the emission.^® A well-established strategy to quantify the amount of spatial fluctu¬ 
ations is the anisotropy APS. Traditionally employed for the study of the CMB [254], 
the technique consists in decomposing a 2-dimensional map /(n) in spherical harmonics 
h£,m(n): /(n) = ^ The APS Q is, then, computed as follows: 


Q 




(24) 


A measurement of the APS of the DGRB provides information on its composition 
that is complementary to the study of its intensity energy spectrum. In Sec. 3.1 we 
summarize the measurement of the DGRB anisotropies performed by the Fermi LAT 
in Ref. [66]. Then, in Sec. 3.2, we describe how such a measurement can be used to 
constrain the DGRB contributors.^^ 


3.1. The Fermi LAT measurement of gamma-ray anisotropies 

The measurement performed by the Fermi LAT Collaboration in Ref. [66] is cur¬ 
rently the only observational data available on the APS of the DGRB. 22 months of 
data are analyzed between 1 and 50 GeV, divided in 4 energy bins. Gamma rays are 
binned into a HEALPix map^^ [402] with N^ide = 512. The count maps are divided by 
the exposure of the instrument in order to obtain gamma-ray intensity maps. This is 
required in order to eliminate any spurious spatial fluctuations due to the non-uniform 
exposure of Fermi LAT. 

Two dehnitions of the APS are used in Ref. [66]. The so-called intensity APS Ci is 
obtained from Eq. (24), while the fluctuation APS comes from the decomposition 

of the relative fluctuations /(n)/(/), where (/) is the all-sky average intensity. The 
two dehnitions are related by = CijlT)^. Note that the huctuation APS is a 

dimensionless quantity, while Ci inherits the units of the intensity map (squared). 

Contrary to what is done in Ref. [9] (see Sec. 2.1), no template htting is employed 
in Ref. [66] to isolate the DGRB. Nevertheless, a mask is applied, screening the regions 
in the sky where the emission is dominated by the diffuse Galactic foreground or by the 
resolved point sources. The mask covers the strip with |6| < 30° around the Galactic 
plane and a 2°-radius circle around each source in the lEGL catalog. The contamination 
of the Galactic foreground is not completely removed by the use of the mask, but the 
residual Galactic emission only induces large-scale features that contribute to the APS 
at small multipoles. That is why only multipoles larger than 105 are considered in 
Ref. [66]. 

The mask may alter the shape and normalization of the APS. Following Ref. [403], 
the APS is corrected for the effect of the mask simply by dividing Ci by the fraction of 


^®Moreover, as we will see in Sec. 3.1, the measurement of the DGRB anisotropies performed in 
Ref. [66] does not rely on the template fitting used in Ref. [9]. 

^®Other observables have been used to quantify the anisotropies in the DGRB: Refs. [107, 398, 399[ 
consider the 2-point correlation function in real space, while Ref. [400] relies on the nearest-neighbor 
statistics. Ref. [401] compares the number of “isolated” gamma-ray events with the “empty regions” in 
the sky. In this section, we focus only on the APS since it is the most commonly used technique. 
http://healpix.jpl. nasa.gov/index.shtml 
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Figure 13: Intensity APS of the data before (purple crosses) and after (red boxes) the removal of the 
diffuse Galactic foreground. The signal region is defined between I = 155 and 504. The large increase 
in the intensity APS with no foreground subtraction (purple crosses) at ^ < 155 is likely attributable to 
contaminations from the Galactic foreground emission. The right panel is an expanded version of the 
left panel and it focuses on the high-multiple angular power. Taken from Ref. [66]. 


unmasked sky, /sky Therefore, the APS estimator considered in Ref. [66] is 


C raw If 

^ ^ £ //sky - Gn 

* ^^^beam^2 ’ 


(25) 


where C'™'" is the intensity APS computed directly from the masked intensity maps. Cn 
is the so-called photon noise, i.e. Cn = (/)^47r/sky/Alj,, where is the total number of 
events detected outside the mask. Finally, the factor corrects for the smearing 

induced by the Point Spread Function (PSF) of the telescope (see Ref. [66] for a dehnition 
of The effect of the PSF becomes too extreme for angular scales beyond (. = 504 

so that, in Ref. [66], no data are considered for multipoles larger than (. = 504. 

Fig. 13 shows the APS for gamma rays between 1.99 and 5.0 GeV. The data points 
show the weighted average of the APS inside bins in multipoles with A£ = 50. The purple 
crosses correspond to the APS estimator Ci of Eq. (25) derived from the intensity maps. 
On the other hand, the empty red boxes indicate the APS measured from the residual 
maps obtained after the subtraction of a model for the diffuse Galactic foreground. The 
error bars are computed following Ref. [404].^^ 

The measured APS is different from zero in the signal region, i.e. 155 < i < 504. 
The (unbinned) APS is htted with a power law (£/155)” in order to determine any 
dependence of the APS on i. R is found that the n = 0 case is compatible with the 
best ht at 95% confidence level for all energy bins. This means that the measured 
APS is compatible with being Poissonian, i.e. independent of i. The significance of the 
detection is 6.5 (7.2) between 1.04 and 1.99 GeV (between 1.99 and 5.0 GeV), while it 
decreases to 4.1 and 2.7 in the two remaining bins at higher energies. 

The fluctuation APS of a population of sources depends on the energy only if their 
spatial clustering is itself energy-dependent or if they are characterized by signihcantly 
different energy spectra. On the other hand, their intensity APS scales with the energy 
as (7)^. When computing the APS of an emission that is the sum of different components 


^^More recently, Ref. [405] presented an alternative method to estimate the error on Ct, developed 
specifically for scenarios like the DGRB, where the emission is affected by limited statistics. 


43 












(as in the case of the DGRB), it may be interesting to study how the APS of the total 
emission depends on the spectra of the individual components By construction, 
the intensity APS is an addictive quantity:^^ 

( 26 ) 

while the fluctuation APS follows the following summation rule: 

^fluct,total _ \ ^ ^fiuct,i /07^ 

— 2-^ Ji ) 

where fi is the fraction of the emission associated with the i-th contribution, with respect 
to the total, i.e. fi = 

Assuming that the components have an energy-independent any energy de¬ 
pendence in arise from the /j-factors in Eq. (27). Indeed, Ref. [359] 

proves that detecting an energy modulation in the fluctuation APS of the DGRB may 
indicate that the emission is the sum of more components (see also Refs. [407, 408]). In 
that case, the behavior of the intensity APS as a function of energy would follow the 
energy spectrum of the dominant component. Thus, studying if (and how) fluctuation 
and intensity APS depend on the energy may be crucial to unravel the composition of 
the DGRB. 

Fig. 14 shows the Poissonian Gp measured in Ref. [66] in the 4 energy bins^^. The 
left panel proves that the fluctuation APS does not depend on the energy.^^ On the 
other hand, the right panel shows that Gp decreases with energy as Gp oc ^-(4.79±o.i3)_ 
This result suggests that the Fermi LAT APS is produced by one single population of 
unclustered sources with an energy spectrum proportional to oc In the next 

section we will see that unresolved blazars fits this description. 

After the publication of the Fermi LAT APS measurement, Ref. [409] pointed out 
that the Fermi LAT Gollaboration used 22 months of data to measure the APS, but 
masked only the contribution of the point-like sources in the IFGL (relative to an ex¬ 
posure of only 11 months). Thus, the emission considered in Ref. [66] will probably be 
contaminated by the contribution of sources that, being not bright enough to be included 
in the IFGL, could have been detected in the larger dataset of Ref. [66]. 

3.2. Deducing the nature of the Diffuse Gamma-Ray Background from its anisotropies 

We start this section by summarizing the technique used to estimate the APS of 
a generic class of sources. The discussion will, then, focus on the gamma-ray emitters 
considered in Secs. 2.2 and 2.3. We will finally show how the comparison of the model 
predictions to the APS signal detected by the Fermi LAT in Ref. [66] can be used to 
constrain the contribution of different populations to the DGRB. 

As done at the beginning of Sec. 2.2, the formalism proposed here assumes that the 
sources are characterized by a generic parameter Y. However, Eq. (1) only determines 


Eqs. (26) and (27) we are neglecting, for simplicity, possible cross-correlation terms between 
different components. Under the hypothesis that gamma-ray sonrces trace the same LSS of the Universe, 
these cross-correlation terms can contribnte significantly to the total APS and, thns, should be 

taken into account [52, 145, 339, 73, 406, 74]. 

^^Since, as commented before, the APS is compatible with being constant in multipole, the whole 
APS can be completely characterized by just one number, that we refer to as Dp. 

However, the large error bars and the fact that only 4 energy bins are available do not allow to 
exclude large-scale modulations or very localized peaks. 
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Figure 14: Anisotropy energy spectra of the data (purple crosses) and of the Galactic-foreground- 
cleaned data (red boxes). Left: Fluctuation anisotropy energy spectrum. Right: Differential intensity 
anisotropy energy spectrum. Taken from Ref. [66]. 


the all-sky average gamma-ray emission associated with population X, referred to here 
as (lx) to underline that it is obtained by integrating over all the possible pointings in 
the sky. The emission /x(n) from direction n can be written as follows [51, 107, 406]: 

Ixitti) = J dxgx{x,i^)Wx{x), (28) 

where x = xi^) is the comoving distance relative to redshift z. Wx{x) is the so-called 
window function and it gathers all the quantities that, in the definition of Ix{n), do 
not depend on the direction of observation. In particular, Wx{x) depend on the 
energy at which Ix is computed. The factor gx{X:^) is called the “source field” and it 
describes how Ix changes from point to point in the sky. It encodes the dependence on 
the abundance and distribution of the sources. The averaged source field, (gx), can be 
used to write (Ix) as f dx(gx(x)}^x(x)- The average source field depends only on the 
abundance of sources and it can be written as;^® 

J min 

The coefficients of the fluctuation APS of population X can be computed 

decomposing the relative intensity fluctuations /x(n)/(/x) in spherical harmonics, as 
follows: 

= (7^/ ^^n/x(n)y^,,(n), (30) 

where dXl^ indicates the angular integration. Now, from Eq. (28) 

j dUn j dxgx{X:T^)Wx{x)Ylm{'^) ( 31 ) 

j dQn j dxfx{x,'n){gx{x))Wx{x), 

where fxiX:'’Y) = gxiX:'’^)/{gx)- The fluctuation APS is obtained from the average of 


fluct,X 


^®Compared to Eq. (1), we are neglecting the dependence on the spectral shape F. 
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with the same multipole i, see Eq. (24). It can be shown [51, 406] that it is 
equivalent to: 




% {9x{x))‘^w\{x) Px (k = -,x) , 

\ X J 


where Pxik,x) is the 3-dimensional power spectrum of the field /x(X)n): 


(/x(x,k)/x(x,k')) = (27r)^(i^(k-F k')Px(A:,x), 


(32) 


(33) 


and fx indicates the Fourier transform of the field fx- Eq. (32) makes use of the so- 
called Limber approximation [410, 411], valid for £ 3> 1. This is indeed the regime of 
interest here since the Fermi LAT APS measurement is robust only for ^ > 155. 

The hypothesis of working with a collection of unresolved sources allows the decom¬ 
position of the 3-dimensional power spectrum Px into a 1-halo term and a 2-halo term, 
Pih and P 2 h, respectively. This helps in the interpretation of Px and it allows to express 
it easily in the context of the Halo Model. The 1-halo term accounts for the correlation 
between two points located within the same source, while P 2 h describes the correlation 
between points that reside in different objects and, thus, it depends on the spatial clus¬ 
tering of the sources. Within the Halo Model formalism, 1-halo and 2-halo terms can be 
written, respectively, as follows [51, 52, 406]: 


Pih{k,z) 


/■Lnaxf^) 

Ir „ 


/ u{k\Y,z)\ 

V iax) J 


(34) 


and 


P 2 h{k,z) = 


'LninD “kdP 


{ax) 


pi“(fc,z). (35) 


The factor u{k\Y, z) is the Fourier transform of the radial brightness profile of a source 
characterized by parameter T at a distance z. Astrophysical sources are normally con¬ 
sidered to be point-like, i.e. intrinsically smaller than the PSF of the telescope at any 
energy. This implies that u is proportional to the source gamma-ray flux S', without 
any dependence on k. In that case, the 1-halo power spectrum becomes Poissonian and 
it depends only on the abundance of sources. In fact, taking Y = S, Eq. (34) can be 
re-written as follows: 


P Poissonian 
Ih 


f 





(36) 


On the other hand, the 2-halo term in Eq. (35) is obtained under the hypothesis 
that the fluctuations in the source distribution traces the fluctuations in the matter 
field, except for a bias factor bxiY^z). That is the reason to include the linear power 
spectrum of matter fluctuations in Eq. (35). Thus, even for point-like sources, the 
2-halo term inherits a dependence on k. 

The balance between the 1-halo and 2-halo terms determines the shape of Pxik,x) 
and of Ci through Eq. (32). For a population of point-like sources that are very bright 
but scarce in number, p^oissoman ^g^^hy dominates and it is difficult (if not impossible) 
to use the APS to extract information on their clustering. On the other hand, if the 
point-like emitters are very abundant, p^issoman jg gj^ialler, allowing for some sensitivity 
to the 2-halo term. On the other hand, if the sources appear extended, the 1-halo is no 
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longer Poissonian and it is suppressed above a certain scale associated with the typical 
size of the sources [51, 107, 406]. 

Following the formalism presented above. Refs. [52, 107] compute the APS of unre¬ 
solved blazars. Blazars are characterized in Refs. [52, 107] by their gamma-ray luminosity 
(y = L^) which is assumed to correlate with the X-ray luminosity Lx, see Sec. 2.2.1. 
The free parameters in the model are fitted to reproduce the abundance of sources de¬ 
tected by EGRET. Ref. [107] finds that the APS of unresolved blazars is within reach 
of the Fermi LAT and that it is dominated by the Poissonian 1-halo term for multipoles 
larger than a few. Similar results are obtained in Ref. [52]. 

Ref. [67] revises the predictions for the APS of unresolved blazars by modeling the 
differential source count distribution dN/dS as done in Ref. [18]. Eq. (36) is used to 
derive the intensity APS from dN/dS: 


Cp 


10 


^thT n AT 


(37) 


where ^thr is the flux sensitivity threshold for point-like sources. dN/dS is assumed 
to be a broken power law, tuned to reproduce the measured abundance of blazars. 
Interestingly, the best-fit dN/dS corresponds to a model in which unresolved blazars 
exhibit an APS that is consistent with the value measured by the Fermi LAT in Ref. [66]. 
This suggests that blazars alone are able to explain the whole APS signal. The same 
best-fit model predicts that unresolved blazars can account only for ~ 30% of the DGRB 
intensity reported by the Fermi LAT in Ref. [8] between 1 and 10 GeV. 

The authors of Ref. [68] employ a more detailed model of the blazar population, 
based on the correlation between and Lx and on a parametrization of their SED (see 
Ref. [19] and Sec. 2.2.1). They confirm the existence of a scenario in which blazars fit at 
the same time the Fermi LAT dN/dS and the measured APS. In this case, unresolved 
AGNs account for at most 4.3% of the DGRB intensity in Ref. [8] (in any energy bin). 
The constraints obtained in Refs. [67, 68] on the contribution of unresolved blazars to the 
DGRB intensity showcase how informative and complementary the study of gamma-ray 
anisotropies can be for the reconstruction of the composition of the DGRB. 

The most recent estimate of the APS of blazars can be found in Ref. [63], where the 
authors consider three distinct AGN subclasses: i) ESRQs, ii) HSP BL Lacs and Hi) a 
combination of LSP and ISP BL Lacs (see Sec. 2.2.1). Their Poissonian intensity APS 
is computed separately, according to Eq. (37). The dN/dS is taken from Ref. [22] for 
ESRQs and from Ref. [24] for the two classes of BL Lacs. They conclude that, as found 
previously, unresolved blazars can explain the whole APS signal, see Eig. 15. HSP BL 
Lacs are responsible for the largest fraction of the measured intensity APS: between 1 
and 2 GeV, they account for 34.5^g4% of the total APS signal and the fraction increases 
to 105]'i3Q% above 10 GeV. 

Ref. [63] also provides the first estimate of APS associated to MAGNs. It is computed 
from Eq. (37), assuming that the 2-halo term can be neglected. The properties of 
MAGNs are inferred from a modeling of the sources at radio frequencies, via the L^ — 
Lr,core relation discussed in Sec. 2.2.2. Unresolved MAGNs are found to contribute to 
approximately 10% of the Fermi LAT APS (6.1% between 1 and 2 GeV and 16.7% above 
10 GeV). 

The APS of unresolved SEGs is computed in Ref. [145]. SEGs are described assuming 
that their luminosity is proportional to the product of their SER and of gas mass. The 
model is tuned to reproduce the properties of the MW. A power law with an index of 
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Figure 15: The angular power C'p(£') in units of E'^'^Cp{AE)~'^ for MAGNs (red long-dashed points), 
a class of sources combining BL Lacs LSPs and BL Lacs ISPs (blue short-dashed), HSP BL Lacs (green 
dotted), FSRQs (yellow dot-dashed) and the total (violet solid) of all radio-loud AGNs. The APS 
measurement by the Fermi LAT Collaboration is also shown (black solid points). Taken from Ref. [63]. 


2.7 is assumed in Ref. [145] as an universal energy spectrum (see Sec. 2.2.3). The APS is 
computed from Eqs. (34) and (35) assuming a bias factor of 1.11, independent of redshift 
and luminosity [412] (see also Ref. [399]). Contrary to the case of blazars or MAGNs, the 
APS of SFGs is not dominated by the Poissonian 1-halo term, at least below multipoles 
of few hundreds. As commented before, this is expected from a population of dim but 
very abundant sources. Thus, the signal may be used to constrain the clustering of 
SFGs. Unfortunately, the signal is overall too faint to contribute significantly to the 
APS detected by the Fermi LAT. 

Refs. [37] and [38] performed Monte Carlo simulations of the APS signal expected 
from unresolved MSPs. Sources are modeled combining radio and gamma-ray data (see 
Sec. 2.2.4). The results of the two references differ from each other, due to the different 
models employed for the description of the MSPs: Ref. [38] assumes that the APS of 
MSPs is Poissonian and it finds that MSPs are responsible for no more than 1% of the 
APS measured by the Fermi LAT, while a larger fraction is allowed in the reference 
model considered by Ref. [37]. 

Among the other astrophysical components considered in Sec. 2.2.5, we mention that 
the clustering of Type la supernovae is considered in Ref. [413]: the authors find that 
unresolved supernovae can exhibit a moderately large APS but the emission peaks in 
the MeV energy range. Refs. [414, 209, 415, 202, 416, 107, 39] study the case of galaxy 
clusters. In particular, the model in Ref. [39] is tested against radio data and it proves 
that the intensity APS associated to galaxy clusters is expected to be two orders of 
magnitude lower than the APS measurement by the Fermi LAT. 

Turning to the case of gamma-ray emission induced by DM, the study of its angular 
anisotropies has also been traditionally considered as a very powerful strategy to sin- 
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Figure 16: Total intensity APS of the gamma-ray emission from DM annihilation (color lines) or decay 
(black line) in extragalactic and galactic (sub)halos. The blue and red lines correspond to the LOW 
and HIGH subhalo boosts, respectively, so that the filled gray area between them corresponds to the 
uncertainty due to the subhalo boost, for a fixed value of Mmin- The red (blue) shaded area around 
the red (blue) solid line indicates the uncertainty in changing the value of Mmin from 10^^ to 1 Mq/H, 
for the LOW (HIGH) case. The red-shaded area is very thin and difficult to see. The solid black line 
shows the prediction for a decaying DM candidate. The observational data points with error bars refer 
to the measurement of the APS between 2 and 5 GeV as given in Ref. [66]. The predictions refer to a 
DM particle with a mass of 200 GeV, an annihilation cross section of 3 x 10“^®cm^s“^ and annihilations 
entirely into b quarks (for annihilating DM) and to a mass of 2 TeV, a decay lifetime of 2 x lO^^s and 
decaying entirely into b quarks (for decaying DM). Taken from Ref. [62]. 


gle out the DM component of the DGRB [417]. Refs. [51, 52] adopt the Halo Model 
approach to determine the 1-halo and 2-halo terms from Eqs. (34) and (35), modeling 
the fluctuations in the gamma-ray emission produced by DM in extragalactic halos and 
subhalos. In Ref. [51] the authors neglect the contribution of subhalos, which is included 
in Ref. [52] assuming that the number of subhalos hosted by a halo with mass M scales 
with M. Their results suggest that, depending on the value of Mmin, tbe DM-induced 
fluctuation APS can be within reach of the Fermi LAT and that it can be larger than the 
APS expected for unresolved blazars.^^ Similar results are obtained also in Ref. [356], 
where only the 2-halo term is considered and modeled from Af-body simulation data. 

The analytic formalism of Refs. [51, 52] is extended to the case of Galactic DM in 
Ref. [54]. The DM halo of the MW is responsible only for large-scale anisotropies. On 


^^The formalism was extended to include p-wave annihilations in Ref. [257]. 
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the other hand, the emission induced by Galactic subhalos can exhibit an APS that 
is even larger than the one from extragalactic DM structures. Results are tested in 
Ref. [54] against different parametrizations of the Galactic subhalo population.^® 

Galactic and extragalactic DM signals are considered at the same time in Ref. [339]. 
For multipoles larger than few tens, the APS of the extragalactic component is dominated 
by the 1-halo term, which depends mainly on the amount of subhalos. In the fiducial 
model of Ref. [339], their luminosity is described by a power-law extrapolation of the 
ca(-2, M) relation. As discussed in Sec. 2.3.1, this probably leads to an overestimation of 
the DM signal, which is, anyway, much smaller than the measured intensity APS. The 
latter can be used to derive upper limits on the annihilation cross section, excluding the 
region in the {m^, (crv)) plane associated with too large anisotropies.®® The upper limits 
obtained in Ref. [339] exclude cross sections larger than 10“^®cm®s“^ (5 x 10“^^cm®s“^) 
for a DM mass of 10 GeV (1 TeV), see the solid red line in Fig. 17. Similar results are 
obtained in Ref. [75]. 

Refs. [53, 55, 359] follows an alternative approach for the study of the anisotropies 
induced by Galactic DM subhalos. Instead of computing the APS analytically like in 
Refs. [51, 52, 54, 356, 339], Refs. [53, 359] simulate sky-maps of the gamma-ray emission 
expected from mock realizations of DM subhalos in the MW. The APS is, then, computed 
by means of the HEALPix package (see also Ref. [419]). In Ref. [366] the results of 
Ref. [53] on the APS of Galactic subhalos are combined with the predictions in Ref. [356] 
for the anisotropies of extragalactic DM to determine the sensitivity of Fermi LAT to 
the detection of a DM component in the DGRB through its APS. 

Similarly, Ref. [55] employs a hybrid approach (inspired by Ref. [274]) in which Galac¬ 
tic DM subhalos are simulated only inside a sphere of radius r^ax centered on the ob¬ 
server. The value of rmax is related to the distance beyond which subhalos become 
point-like. It is assumed that subhalos located further than rmax cumulatively generate 
a smooth emission. This is equivalent to assume that APS generated by subhalos beyond 
^'max is dominated by the 2-halo term. 

Semi-analytic hybrid methods are used also in Refs. [56, 62] to compute the anisotropies 
in the emission of extragalactic DM structures. In this case, the distribution and prop¬ 
erties of extragalactic (sub)halos with a mass larger than the mass resolution of the 
Millennium-II A^-body simulation [343] are taken directly from the halo catalogs of the 
simulation. Mock sky-maps are generated by replicating the Millennium-II simulation 
box until it covers the region within z ~ 2.^® DM halos less massive than the resolution 
of Millennium-II are included assuming that they share the same clustering of those 
immediately above the mass resolution [56, 343]. Subhalos of extragalactic DM clumps 
are also accounted for considering multiple scenarios for their abundance and internal 
properties. Ref. [62] completes the prediction by modeling also the Galactic signal. The 


^®Ref. [54] and the majority of the references mentioned in this section do not consider the effect of 
baryons on the clustering of DM and, thus, on its anisotropy pattern. One exception is Ref. [418], which 
computes the APS expected by the so-called DM “mini-spikes”, i.e. DM overdensities induced by the 
presence of Intermediate-Mass Black Holes. The authors find that this scenario can lead to a significant 
increase in the amplitude of the DM-induced APS. 

^®We note that, in general, the intensity APS associated with DM scales quadratically on {av). The 
study of anisotropies is, therefore, quite sensitive to a DM signal. Whether the upper limits on {gv) 
derived from the Fermi LAT APS are more constraining than those from the DGRB intensity energy 
spectrum, will depend on the amplitude of the intrinsic fluctuations in the DM-induced gamma-ray 
emission. See also Fig. 17. 

"^^Beyond this distance, at the energies considered, the gamma-ray flux is almost entirely attenuated 
by the interaction with the EBL. 
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all-sky DM maps produced in Ref. [62] represent complete and accurate templates of 
the total DM-induced gamma-ray emission. 
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Figure 17: Upper limits obtained by considering the DGRB APS measured in Ref. [66]. The regions 
above the solid red and blue lines are excluded because the cumulative DM-induced emission would 
overproduce the APS of the DGRB. The solid red line is taken from Ref. [339] while the solid blue 
one is from Ref. [420]. The dashed gray line show the upper limits on (av) obtained by requiring that 
the DM-induced emission does not overproduce the DGRB emission measured in Ref. [9]. The limit is 
obtained in Ref. [25] by modeling the astrophysical DGRB contributors. The blue region indicates the 
region already excluded by the observation of the Segue 1 dwarf Spheroidal galaxy performed by the 
MAGIG telescope [368]. The dark gray region is excluded by the analysis performed by the H.E.S.S. 
telescope in Ref. [369] on the so-called “Galactic Center halo” region (assuming an Einasto DM density 
profile). The light gray region indicates the DM candidates not compatible with the combined analysis 
of Fermi LAT data from 15 dwarf Spheroidal galaxies [370]. The dash-dotted line marks the thermal 
annihilation cross section. 



Ando and Komatsu (2013) MAGIC Segue 1 

Gomez-Vargas et al. (2014) H.E.S.S. Galactic Center halo 

Ajello et al. (2015) i-r-Tij Fermi LAT dwarf Spheroidals 


As in Sec. 2.3.1, the main sources of uncertainties in the DM-induced APS are the 
value of Mmin and the subhalo boost factor. Extreme scenarios are identified in Ref. [62] 
for both Minin and the subhalo boost factor, bracketing their theoretical uncertainties. 
The effect of their variability on the intensity APS sums up to approximately two orders 
of magnitude, as it can be seen in Fig. 16. The red and blue lines show the intensity 
APS for two different prescriptions of the subhalo boost, labeled “LOW” and “HIGH”. 
The blue line is obtained following the prescription of Ref. [276], which employs power- 
law extrapolations for low-mass halos. On the contrary, the red line makes use of the 
model by Ref. [327], extended in Ref. [329]. The blue-shaded region shows the effect of 
changing on the HIGH scenario.The black data points indicate the Fermi LAT 


Changing Mmin has no effect on the red line (LOW case), see Ref. [62] for details. 
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APS measurement between 2 and 5 GeV. As in Ref. [339], the DM-induced emission is 
found to contribute only marginally to the measured APS. 

In Ref. [420] the predictions of the HIGH case (blue line) are, then, used to derive 
upper limits on the DM annihilation cross section, shown as a solid blue line in Fig. 17. 
The upper limit reaches values as low as 3 x 10“^®cm^s“^ for m-^ = 50 GeV, while going 
to 10“^^cm^s“^ if the mass approaches the TeV scale. Fig. 17 also summarizes other 
constraints on {av) available in the literature. In particular, we plot the limit from the 
analysis of Ref. [25] of the DGRB intensity energy spectrum (dashed gray lines). The 
shaded regions indicate the areas already excluded by other indirect searches of DM, 
i.e. the observation of Segue 1 with MAGIG (blue), of the Galactic Center halo region 
with H.E.S.S. (dark gray) and the Fermi LAT combined analysis of dwarf Spheroidals 
(light gray). Fig. 17 proves that the study of the anisotropies of the DGRB can produce 
competitive upper limits on {av), but probably not as strong as those induced by the 
DGRB energy spectrum. 

Ref. [62] also estimates the APS associated with a decaying DM candidate. In this 
case, as found also by Ref. [57], the predictions are subject to less theoretical uncer¬ 
tainties than for an annihilating DM candidate. In fact, the signal is less affected by 
the value of Mmin and there is no subhalo boost (see also Sec. 2.3.2). Yet, in decaying 
DM scenarios, DM halos yield a more extended emission. This is particularly true for 
Galactic subhalos which are still close enough not to be point-like. Thus, the APS is 
expected to decreases rapidly at high multipoles being, therefore, hard to detect. See, 
e.g., the black line in Fig. 16. 

4. The photon count distribution 

Another powerful statistic tool to constrain the nature of the DGRB is provided 
by the photon count Probability Distribution Function (PDF). This technique can be 
used when the emission is represent by a pixelated sky-map. The photon count PDF 
is, then, built from the number of pixels in which k photons are detected. The 
study of the photon count PDF is commonly used in radio and X-ray astronomy for 
the analysis of diffuse emissions, in particular when trying to estimate the contribution 
of faint unresolved sources. It can also be used at gamma-ray frequencies to single out 
different components in the DGRB, even if they are subdominant. Indeed, different 
PDFs are expected for different populations of gamma-ray emitters: bright but rare 
sources generate pixels with a large (or moderately large) number of photons, i.e. Uk 
will be significantly different from zero even at large k. On the other hand, a population 
of faint but numerous sources is normally associated with a Poissonian PDF. Intrinsically 
diffuse emissions also correspond to Poissonian PDFs. 

Ref. [69] measures the PDF directly from 11 months of Fermi LAT data between 1 
and 300 GeV. The data are binned into a HEALPix map with Aside = 32, corresponding 
to a binsize of approximately 0.4°. A region of 30° around the Galactic plane is masked, 
while point sources are not masked. The observed PDE is represented as red data 
points in Eig. 18. A model based on the so-called generating functions is developed to 
interpret the data. If pk is the probability of finding k photons in a certain pixel (i.e. 
Pk = ^^fc/Apixeis, where Apixeis is the total number of pixels considered in the analysis), 
the corresponding generating function P{t) is defined as 

OO 

P{t) = Y^PktK (38) 

k=0 
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Viceversa, the probabilities pk can be derived from P{t) as the coefficients in front of 
each term in the power-law expansion of P{t). The generating function of an emission 
composed by multiple components is the product of the generating functions of the single 
components [421]. In particular, the model considered in Ref. [69] consists of 3 terms: 

• point sources, both resolved and unresolved. The corresponding generating func¬ 
tion is P{t) = — Xm)], where Xm is the average number of sources 

emitting exactly m photons in a pixel. The set of {xm} can be derived from the 
differential source count distribution dN/dS. In Ref. [69], dN/dS is taken to be a 
broken power law [18], leaving its 4 parameters free. 

• the diffuse Galactic foreground, whose emission, in each pixel, is interpreted as a 
collection of sources that emit exactly 1 photon. The number of those sources is 
proportional to the intensity of the foreground in that pixel. Its generating func¬ 
tion P{t, i) depends on the particular pixel i considered and it can be written as 
P{t, i) = exp(x[j;gt—x[j;g), where is the number of photons of foreground emis¬ 
sion in the i-th pixel. In Ref. [69], the intensity and morphology of the foreground 
are derived from the Fermi LAT data itself. 

• an isotropic component, included to represent the emission associated to faint 
abundant sources. As for the Galactic foreground, this term is modeled as a 
collection of xiso sources emitting exactly 1 photon and its generating function is 
P{t) = exp(xisot — a^iso), for ah the pixels. The quantity Xiso is left free in the fit. 

The model is fitted to the photon count PDF in Fig. 18 in order to determine the 
free parameters, i.e. those characterizing dN/dS and the normalization of the isotropic 
component. The best-fit dN/dS is compatible with what found in Ref. [18]. The analysis 
of the PDF is different and complementary to Ref. [18] but it succeeds in reconstructing 
the properties of a population of non-Poissonian blazar-like sources. By integrating the 
best-fit dN/dS obtained in Ref. [69] below the sensitivity of the Fermi LAT, unresolved 
blazars are found to account for ~ 23% of the DGRB reported in Ref. [8] above 1 GeV. 
We remark that this result is consistent with the independent estimation obtained in 
Ref. [18]. Similar results are also obtained in Ref. [364] using 5 years of simulated Fermi 
LAT data and a pixel size of 0.25°. 

Predictions for the PDF are available not only for blazars: Ref. [36] computed the 
PDF expected from unresolved MSPs finding that their PDF is highly non-Poissonian. 

In the case of DM-induced emission, the shape of the PDF is expected to depend 
on which DM structures are considered. Ref. [422] computes the probability P{F) 
of detecting a certain flux F due to Galactic DM subhalos in a generic pixel. The 
probability depends on the modeling of the subhalo population and, in particular, on 
the value of Mmin- It can be written as follows (see the Appendix of Ref. [422] for a 
detailed derivation): 

P{F) = F-^{eMl^{F{Pi{F)] - 1)]}, (39) 

where F{f{x)} indicates the Fourier transform of /(x) and its inverse. The quantity 
p is the average number of sources expected inside one pixel and Pi{F) is the probability 
of having exactly one source emitting the flux F. The DGRB model in Ref. [422] also 
includes a background component with a Poissonian PDF. The authors prove that there 
exists a region in their parameter space in which Galactic DM subhalos are faint enough 
to escape detection as individual sources after 5 years of Fermi LAT data, but bright 
enough to be detected by studying the photon count PDF. 
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Figure 18: Uk is the number of pixels with k photons. The red dots correspond to the pixel counts 
derived from the data of the Fermi LAT, the error bars are equal to ydifc- The model devised to 
interpret the data considers three classes of sources: i) AGN-like point-like objects (blue dotted line), ii) 
isotropic Poisson contribution (brown dashed line) and Hi) anisotropic Galactic diffuse emission (black 
dash-dotted line). Note that the photon count PDF for the total (solid black line) on this plot is not 
the sum of the components, but the corresponding generating function of the PDF is the product of 
the generating functions of the three contributions, ^parameters is the number of parameters used in 
the fit. qps, Qiao and Qgai are the relative contributions of the point sources, of the isotropic and of the 
Galactic foreground, respectively. Appoints is the number of points in the a;-axis employed in the fit and 
the log-likelihood of the best-fit is indicated in the legend. Taken from Ref. [69]. 


Ref. [423] also studies the possibility of separating the emission associated with Galac¬ 
tic DM subhalos from the diffuse Galactic foreground and from an isotropic background 
(both characterized by a Poissonian PDF). The authors choose some benchmark mod¬ 
els for the description of the DM-induced emission and of the Poissonian backgrounds. 
After building the PDFs of the different components, they use them to generate mock 
sky-maps of the expected gamma-ray emission. This mimics a real observation and it 
allows to estimate the sensitivity of a PDF analysis to the DM signal. The mock photon 
count PDF is compared to model predictions in order to determine the free parameters 
in the model, as done in Ref. [69]. This relies on a a priori knowledge of the shape of the 
PDF for the different components in the model. Ref. [423] proves that, in an idealized 
case without background, the reconstructed intensity of the DM-induced emission de¬ 
pends significantly on the assumed shape of the PDF. Employing the wrong PDF would 
lead to biased reconstruction of the intensity of the DM component. However, when the 
Poissonian backgrounds are included. Ref. [423] finds that an unbiased reconstruction 
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(within statistical uncertainties) can be obtained independently of the shape of the PDF. 

Ref. [424] considered the gamma-ray emission of extragalactic DM halos and sub¬ 
halos, instead. The authors derive Pi{F) analytically, relying on the halo model, and 
they show how its shape changes when using three different prescriptions for the subhalo 
boost. Then, they compute P{F) by means of the central limit theorem below a refer¬ 
ence flux T*, and through MC simulations above that. The resulting flux distribution 
P{E) exhibits two regimes: below 5 x 10“^^cm“^s“^GeV“^sr“^ it follows a Gaussian 
distribution, while it is a power law for larger fluxes.The latter is the case when a few 
bright DM structures dominate the flux distribution. Ref. [424] also developed a model 
for the P{F) of the astrophysical components of the DGRB: for a fiducial subhalo boost 
model inspired by Ref. [280], the authors of Ref. [424] show that a measurement of {pk} 
after 5 years of Fermi LAT data can lead to a detection of a DM signal, provided that 
the annihilation cross section is, at least, twice the thermal value, for a DM mass of 85 
GeV and annihilations into b quarks. 

5. The cross-correlation with independent probes 

The most recent development in our understanding of the composition of the DGRB 
has focused on the study of its cross-correlation with other observables. For instance, the 
fraction of the DGRB that originates from extragalactic objects (whether astrophysical 
sources or DM halos) traces the LSS of the Universe, up to a maximal redshift that 
depends on the EBL attenuation. Thus, a certain level of cross-correlation is expected 
with any LSS tracer, e.g. the distribution of resolved galaxies [399, 73, 74, 70] or the 
gravitational leasing effect of cosmic shear [72, 71, 75]. These are new and indepen¬ 
dent observables that can complement the information inferred from the DGRB energy 
spectrum (see Sec. 2) or from its auto-correlation APS (see Sec. 3). Also, note that 
the gamma-ray sky-maps analyzed in Ref. [66] to measure the auto-correlation APS are 
noise dominated. Thus, even if it was still possible to report a significant auto-correlation 
signal by subtracting the photon noise, one may expect the cross-correlation with other 
signal-dominated quantities to be, in principle, very informative. 

In the following sections, we present the formalism proposed to predict the cross¬ 
correlation of the DGRB with LSS tracers. We also summarize the data currently 
available. Sec. 5.1 focuses on the cross-correlation with galaxy catalogs, while in Sec. 5.2 
we discuss the case of the cosmic shear. Sec. 5.3 presents the results of the cross¬ 
correlation of the DGRB with other observables. 

5.1. The cross-correlation with galaxy catalogs 

Ref. [399] was the first work that measured the 2-point correlation function of the 
Fermi LAT DGRB with 4 galaxy surveys, namely: i) the Sloan Digital Sky Survey 
(SDSS) Data Release 6 of optically-selected quasars from Ref. [425], ii) the IR galaxies 
of the 2 Micron All-Sky Survey (2MASS) Extended Source Catalog from Ref. [426], Hi) 
the radio sources from the NRAO VLA Sky Survey (NVSS) [427] and iv) the luminous 
red galaxies in SDSS Data Release 8 from Ref. [428]. The authors of Ref. [399] analyzed 
21 months of Fermi LAT data but no significant cross-correlation was observed. More 
recently. Ref. [70] have updated the analysis using 60 months of data and exploring 
also the cross-correlation with the main galaxy sample of SDSS Data Release 8 from 
Ref. [429]. The region with |6| < 30° around the Galactic plane was masked out, as well 


■^^The gamma-ray flux is computed at 1 GeV. 
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as the 1-degree region around the point sources in the 3FGL catalog. Both the Large 
and Small Magellanic Clouds were also left out of the analysis, together with the Fermi 
Bubbles and of the so-called Loop-I. A model for the diffuse Galactic foreground was 
subtracted from the sky maps. The residuals of the gamma-ray maps and the distribution 
of galaxies in the catalogs were re-binned into HEALPix maps with Aside = 512. 

The authors computed both the cross-correlation APS and the 2-point correlation 
functions^^ using the PolSpice package^^. The signal region was defined between 0.1 
and 100 degrees for the 2-point correlation function and between multipoles of 10 and 
1000 for the cross-correlation APS. Three different energy thresholds were considered 
for the gamma rays, including all events above i) 500 MeV, ii) 1 GeV or in) 10 GeV. 
The authors in Ref. [70] also validated their results against changes in the mask and in 
the model adopted for the diffuse Galactic foreground. They also tested their analysis 
pipeline on a simulated sky map with no signal. 

The 2-point correlation functions are shown in Fig. 19 for the 5 different catalogs: 
a cross-correlation signal is evident up to few degrees for all the catalogs apart from 
the SDSS luminous red galaxies (rightmost central panel). The significance of these 
detections is 4.5a for the SDSS quasars, 3.6 cj for the 2MASS catalog, 3a for the SDSS 
main galaxies and as large as lOci in the case of NVSS. These numbers refer to the 
energy threshold that maximizes the detection significance, i.e. 500 MeV for the case 
of the SDSS quasars and 1 GeV otherwise. These are also the thresholds considered in 
the different panels of Fig. 19. The case of the cross-correlation with NVSS deserves 
some additional comments: the authors of Ref. [70] noticed that, for that catalog, the 
cross-correlation signal exhibits an angular extension that is consistent with the PSF of 
Fermi LAT. In particular, it decreases as the energy threshold increases. This suggests 
a different origin for the NVSS signal with respect to the one observed with the other 
catalogs. They also noted that a 1-halo component to the signal (see Sec. 3.2) would 
manifest itself as a Dirac delta at 0 = 0 degrees, smeared up to the angular size of 
the Fermi LAT PSF. Such a 1-halo term would be present, for example, if some of the 
galaxies in NVSS emitted also in the gamma-ray band as well. Indeed, NVSS galaxies 
are standard candidates to be gamma-ray emitters and this catalog is routinely searched 
for counterparts of gamma-ray sources [115, 2]. With all these considerations in mind, 
the authors of Ref. [70] concluded that the 2-point correlation function with NVSS is 
probably contaminated by a 1-halo term and it does trace the LSS. 

In Ref. [70], the measured correlation functions are compared with the predictions ob¬ 
tained if the DGRB were contributed entirely by only one class of astrophysical sources. 
These theoretical predictions are included in Fig. 19 as colored lines. In particular, 
the dashed red line stands for FSRQs, modeled as in the luminosity-dependent density 
evolution scheme of Ref. [22], while the solid black is for BL Lacs, following the results 
of Ref. [23], see Sec. 2.2.1. The dot-dashed lines denote the case of SFGs, considering 
two different models: i) the gamma-ray emission of SFGs is assumed to follow the cos¬ 
mic SFR, as in Ref. [31], and ii) the gamma-ray luminosity is related to the IR one 
through the L.y — Ljr relation obtained in Ref. [33] (see Sec. 2.2.3). The two scenarios 
are represented by the blue and green dot-dashed lines, respectively. Finally, we note 
that MAGNs are not included in the analysis of Ref. [70] as their emission is showed 
to be very similar (and, therefore, degenerate) with the contribution of SFGs. MAGNs, 


note that the two are related by a Fourier transform and they contain the same information. 
However, they probe different scales with different efficiency and, thus, it is interesting to consider both. 
http://www2.iap.fr/users/hivon/software/PolSpice/ 
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Figure 19: Upper left: 2-point cross-correlation fnnction (orange data points) estimated from the SDSS 
Data Release 6 optically-selected quasars and the DGRB obtained from 60 months of Fermi LAT data 
at l&l > 30° and for E > 500 MeV. Errors bars represent the diagonal elements of the covariance matrix. 
Model predictions for different classes of sonrces are represented by continuons cnrves: FSRQs (red 
dashed), BL Lacs (black solid), two models of SFGs (blue and green dot-dashed). All predictions are 
obtained assuming that each source class contributes 100% of the DGRB intensity and they do not 
represent hts to the data. Upper right: same as in the previous panel but for the cross-correlation with 
the 2MASS Source Extended Catalog and for energies l 1 GeV. The other panels show the same as in 
the previous panel but for the NVSS catalog (medium left), luminous red galaxies in the SDSS Data 
Release 8 (medium right) and the main galaxy sample in the SDSS Data Release 8 (bottom). Note the 
different scale in the plots. Taken from Ref. [70]. 
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however, are explicitly included in the follow-up analysis of Ref. [430] . 

The 2-point correlation functions are, then, computed assuming that the fluctuations 
in the gamma-ray maps and in the galaxy distributions from the catalogs both trace the 
LSS matter density fluctuations (provide that the so-called bias factor is considered). 
This allows for the correlation functions to be determined in terms of the non-linear 
power spectrum of matter fluctuations, which can be obtained, e.g., by means of the 
public CAMB code [431] or Halofit routine [432]. 

We note that the cross-correlation expected from the classes of astrophysical sources 
mentioned above (colored lines in Fig. 19) is indeed different from zero at small angles, as 
in the observed signal. Whether the amplitude of the predicted correlation functions is 
in agreement with the data or not depends on the redshift overlap between the gamma- 
ray emitters and the sources in the catalogs. In particular, optically-selected quasars 
and sources in the NVSS catalog have a quite broad redshift distribution, extending to 
2 ; ~3-4. Thus, a large cross-correlation is expected with the emission from unresolved 
SFGs (modeled as in Ref. [31]), whose redshift distribution peaks at 2 ; ~ 2 — 3. On the 
other hand, the luminous red galaxies and the objects in the 2MASS catalog probe the 
local Universe (respectively, from 2 ; ~ 0.8 and 2 : ~ 0.3, and down to the present time) and 
they are characterized by narrower redshift distributions. Both are expected to correlate 
mainly with BL Lacs. Indeed, by considering several galaxy catalogs with different 
redshift distributions, one can effectively probe different redshift ranges, developing a full 
tomographic approach. Specifically, in Ref. [70], the authors build a model of the DGRB 
that includes FSRQs, BL Lacs and SFGs, leaving the normalization of the different 
components free to vary when fitting the cross-correlation data in Fig. 19. Ref. [70] 
shows that including the cross-correlation with SDSS quasars is crucial in deriving a 
lower bound on the contribution of SFGs."^^ These are, indeed, found to be the dominant 
component in the DGRB, with blazars accounting for, at most 10% of the total DGRB 
measured by Fermi LAT in Ref. [9] (at Icr level). Also, depending on which description 
is assumed for the SFGs, the best-ht model to the cross-correlation data in Fig. 19 can 
account for only 70% or 20% of the total DGRB intensity. 

The cross-correlation with galaxy catalogs can also be used to constrain the DM 
component of the DGRB. This possibility was studied in Refs. [73, 74], where the authors 
considered the 2MASS Redshift Survey [433] and the 2MASS Extended Source Catalog 
from Ref. [426]. These two catalogs are chosen against others because they trace the 
matter distribution in the local Universe and, therefore, they are expected to correlate 
with any potential DM-induced gamma-ray signal. The 2MASS Redshift Survey extends 
to 2 ; ~ 0.1, while the 2MASS Extended Source Catalog peaks at 2 : = 0.072 and does not 
contain galaxies beyond 2 ~ 0.4. The Fermi LAT has detected most of the blazars in 
this volume, down to a very low sensitivity. Thus, unresolved blazars are not expected to 
exhibit a large cross-correlation with catalogs of the local Universe.^® This means that 
the cross-correlation APS will be potentially very sensitive to other DGRB contributors, 
emitting mainly at low redshift as, e.g., SEGs or DM. 

In Refs. [73, 74], galaxies are described by means of the so-called Occupation Distri¬ 
bution (HOD) model. This framework postulates that each source is embedded into a 
DM halo of mass M and that the abundance and distribution of galaxies are related to 
the properties of the host DM halos. The HOD model is a phenomenological formalism, 


"^®Note that this lower bound is found only if SFGs are modeled according to Ref. [31]. 

'^®In other words, the window function of unresolved blazars does not overlap significantly with that 
of the 2MASS catalogs considered above. 
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Figure 20: Left: la and 2cr allowed regions for the DM annihilation rate versus DM mass, for different 
gamma-ray production channels and assuming the “LOW” substructure scheme in Ref. [62]. Crosses 
indicate the best-fit models. In the “HIGH” scenario of Ref. [62], regions remain similar in shape but 
they shift downward by a factor of ~ 12 (star symbols). Right: The same but for decaying DM, showing 
the DM particle lifetime as a function of its mass. Taken from Ref. [436]. 


based on results from iV-body simulations and on semi-analytical descriptions of DM 
halos and galaxy formation [434, 201, 435]. Under this formalism. Ref. [73] shows that 
DM dominates the cross-correlation with the 2MASS Redshift Survey and that 5 years 
of Fermi LAT data should be able to distinguish a DM scenario from one with purely a 
strophysical sources. Ref. [74] also computes the upper limit that it would be possible 
to derive on (av) if the cross-correlation were found compatible with an astrophysical 
interpretation. These results, however, largely depend on the model adopted for low- 
mass DM halos and subhalos. In the most optimistic scenario, the cross-correlation will 
be able to exclude thermal cross sections for DM masses up to almost 1 TeV. 

These predictions were tested against actual data in Ref. [436], where the authors 
explained the measured cross-correlation signal with the 2MASS Extended Source Cat¬ 
alog in Ref. [70] in terms of DM. They found that a WIMP DM particle with a mass 
in the range between 10 and 100 GeV (depending on the annihilation channel) and a 
thermal cross section can reproduce both the shape and intensity of the measured cross 
correlation. This is also the case for decaying DM candidates, for a similar range in DM 
mass and a decay lifetime between 5 x 10^® and 5 x 10^^ s. The colored contours in 
Fig. 20 show the the regions in the {my-,{(Tv)) and {m^,T) parameter spaces that are 
compatible with the cross-correlation found with 2MASS. The figure also shows how, 
for an annihilating DM candidate, a different assumption for the subhalo boost can shift 
the preferred contours by more than one order of magnitude. 

In additiona. Ref. [436] used the measured cross-correlation to derive upper limits 
on {(Tv), by requiring that the DM-induced correlation function do not to over-produce 
the data. This allowed Ref. [436] to exclude thermal annihilation cross sections for DM 
masses below 100 GeV (in the case of annihilations into b quarks and a “LOW” subhalo 
boost model inspired by Ref. [62]). This makes the cross-correlation with local galaxy 
catalogs the strongest observable up to date to constrain a potential DM contribution 
to the DGRB, compared to the DGRB energy spectrum reported in Ref. [9] or the 
auto-correlation APS in Ref. [66]. 

Astrophysical and DM-induced emissions are considered at the same time in Ref. [430] : 
the authors define a model of the DGRB that includes FSRQs, BL Lacs, SFGs, MAGNs 


59 










1 


0.5 


0 

6 




in 

u. 


< 


2 


0 

0.6 


1/1 

O 0.4 
u. 


< ' 0 2 
0 
1.5 


0.5 


0 0.2 0.4 0.6 



0 0.2 0.4 0.6 0 2 4 6 0 0.2 0.4 0.6 0 0.5 1 


^ 2 
n 


< 0 




Figure 21: Posterior probability distributions for the normalization of the component to the DGRB 
due to FSRQs (^fsrQs), BL Lacs (^BLLacs), SFGs (Asfgs), MAGNs (AmAGNs) and annihilating DM 
(ArmDM)- These parameters are defined with respect to a fiducial model introduced in Ref. [430]. Panels 
along the diagonal show the marginalized 1-dimensional probability distribution for each parameters. 
All the others indicate the Ict (darker, innermost) and 2a (lighter, outermost) confidence level contours 
in the probability distributions of the different combinations of parameters. The full model contains 
a total of 11 free parameters, but only the 5 mentioned above are shown in this figure. Taken from 
Ref. [430]. 
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Figure 22: 95% confidence limit upper bounds on the DM annihilation rate (cru) as a function of the 

DM mass, for the “LOW” substructures model of Ref. [62] and the reference NVSS-10 7^ 0 fit (see 
Ref. [430] for details). Solid lines refer to the bb annihilation channel: the red line refers to the analysis 
that combines information from all the three energy bins under consideration (i.e. E > 0.5, 1 and 10 
GeV), while the other three lines refer to the analysis performed on a single energy bin (as stated in 
the figure label). The upper dot-dashed blue line refers to the “NS” substructure model, where halos do 
not have substructures and Mmin ~ IO^Mq. The lower dot-dashed black line, instead, representto the 
“HIGH” substructure model, inspired by Ref. [62]. Taken from Ref. ]430]. 


(parametrized according to Refs. [22, 23, 161, 29], respectively) and annihilating DM. 
The galaxy catalogs are described following the HOM formalism. The normalizations of 
the emission of the 4 mentioned astrophysical source classes are left free in the model, 
as well as the DM mass and annihilation cross section. 5 additional parameters are 
included, one for each class of gamma-ray emitters, accounting for possible 1-halo terms 
in the 2-point correlation functions. The model is, then, used to fit the measured cross¬ 
correlation reported in Ref. [70]. As expected, the posterior probability distribution 
function for the intensity of the 1-halo term points towards a value different than zero, 
in the case of the cross-correlation with NVSS. The distributions are compatible with 
zero for the other catalogs. 

Fig. 21 summarizes the probability distributions of the remaining 5 free parameters 
(excluding m^). An examination of this figure makes evident that the interpretation of 
the cross-correlation data is affected by significant degeneracies: there is a mild indica¬ 
tion of a peak in the probability distribution of the normalizations of MAGNs (AmAGNs 
in Fig. 21) and of the DM component (Admj proportional to {crv)), but, otherwise, the 
data effectively just enforce upper limits on the other parameters. The degeneracy is 
particularly visible between AmAGNs and A^m in the fourth panel of the bottom row of 
Fig. 21. Nevertheless, the measured cross correlations are still able to provide stringent 
uppers limit on the DM component, as it can be seen in Fig. 22. The regions above 
the solid lines are excluded as they would over-produce the 2-point correlation functions 
measured above 500 MeV (yellow line), 1 GeV (cyan line), 10 GeV (green). The red 
solid line indicates the upper limit obtained from the combined analysis of the data sets 
for the three mentioned energy thresholds. These results are all obtained by adopting 
the assume a “LOW” subhalo boost. Note that, even for such a moderate description of 
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low-mass DM structures, the inferred DM limits are more stringent than those obtained 
in Ref. [64] by studying the DGRB intensity or those derived from the DGRB auto¬ 
correlation APS in Ref. [420].^^ In the case that a “HIGH” subhalo boost were adopted 
instead, the upper limit on the DM annihilation cross section would improve by roughly 
one order of magnitude (dash-dotted black line in Fig. 22). On the contrary, in an overly 
conservative scenario in which no substructures are present below Mmin = IO^Mq, the 
exclusion worsen by a factor ~ 4 — 5 (blue dash-dotted line). 

As already hinted at in Ref. [70], the authors of Ref. [430] noted that the models that 
provide a good ht to the cross-correlation data fall short of accounting for the intensity 
of the DGRB: above 1 GeV, in the most likely scenario, the model is able to explain only 
~ 30% of the DGRB intensity measured by Fermi LAT in Ref. [9]. Note, though, that 
multiple reasons can be invoked to explain this apparent discrepancy, e.g. uncertainties 
in the modeling of the Galactic diffuse foreground when performing the measurement of 
the DGRB intensity, and/or uncertainties in the model predictions. Another possibility 
is the presence of a Galactic component in the DGRB that, therefore, does not correlate 
with the LSS. In this regard. Galactic DM would be a plausible and interesting candidate 
(see, e.g.. Ref. [64] for further discussion on this issue). 

5.2. The cross-correlation with cosmic shear 

Another tracer of LSS is the gravitational leasing effect of cosmic shear. Due to 
leasing, the light emitted by distant sources is distorted while it propagates towards us 
by the presence of intervening matter. In the weak leasing regime, the effect is very small 
and it is directly related to the distribution of matter at large scale. The signal, referred 
to as cosmic shear, is expected to cross-correlate with the gamma-ray emission since the 
same structures responsible for light bending are also those producing the gamma-ray 
emission, either because they emit light themselves (through DM annihilation or decay) 
or because they host the astrophysical emitters. 

The gravitational distortion can be evaluated on the null-geodesic of the unlensed 
photon. It can be decomposed into the so-called convergence k and shear 7 [437, 438]. 
In the flat-sky approximation (i.e. small distortion angles), k and 7 share the same APS 
and, for convenience, we focus only on the former from now on. The convergence k is 
a direct estimator of the fluctuations in the Newtonian potential of the LSS, integrated 
along the line of sight, k can be estimated via a statistical analysis of the correlations 
in the ellipticities of the images of galaxies. Thanks to Poissons equation, which links 
the gravitational potential to the matter distribution, the intensity of the cosmic shear 
signal /^(n) from a direction n can be written as follows: 

/^(n) = j dx9K{n,x)Wi^ix)- (40) 

The decomposition of 7k( n) resembles the way we expressed the gamma-ray emission in 
Eq. (28): the source field ^^(n, x) indicates directly the distribution of matter and it is 
modulated by the window function. The product of the average source field {gnix)) and 
the window function can be expressed as follows: 

f°° v' — V dN Iv'i 

{gK{x))W^ix) = +z{x)]x J - dx^' 


47 


These limits are sensitive to the way the 1-halo terms are implemented. See Ref. [430] for details. 
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where dNg/dx' is the redshift distribution of the background galaxies, normalized to 1 
over the observed redshift range. 

As done in Sec. 3, the gamma-ray emission associated with a certain population X 
is written in terms of its window function Wx{x) of the source field gxix^'^ 
/x(n) = f dxgx(Xj n)^x(x)- The average source field (gx(x)} depends on the abun¬ 
dance of sources as a function of their y-parameter and of redshift (see Eq. (29)). Sim¬ 
ilarly to Eq. (32), the cross-correlation APS between the gamma-ray emission produced 
by population X and the cosmic shear can be written as [406, 74]: 


The 3-dimensional cross-correlation power spectrum Px,k can be split into the fol¬ 
lowing 1-halo and 2-halo terms: 


Pih(k,z) 
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(44) 

The quantity u{k\M) is the Eourier transform of the radial density profile of a DM halo 
with mass M, while ux{k\Y{M)) is the Eourier transform of the gamma-ray surface 
brightness profile of the source characterized by parameter Y. Note that these equations 
depend on the Y{M) relation that links the T-parameter to the mass of the host DM 
halo. Eor astrophysical sources, Y is usually taken to be the gamma-ray luminosity L.y. 
Ref. [75] determines the L.y(M) empirically for the case of blazars, SEGs and MAGNs, 
making use of correlations between different source properties or of the results of semi- 
analytical models. However, the L^{M) remains very uncertain for all the source classes 
considered in Ref. [75]. Its uncertainty can become an issue when estimating the cross¬ 
correlation APS However, in the case of astrophysical sources, Cf’'^ is dominated 

by the 2-halo term (at least up to multipoles as large as few hundreds), while the L^{M) 
relation mainly affects the 1-halo term.^® The uncertainty on generated by the 

variability of L^(M) can be seen in the left panel of Pig. 23. The different solid lines 
show the expected cross-correlation APS for different classes of astrophysical sources 
(red for blazars, orange for SEGs and pink for MAGNs), while dashed lines indicate 
how the cross-correlation APS changes for two extreme scenarios bracketing our lack of 
knowledge on XI{L^). The uncertainty bands are within a factor of 2 from the solid 
lines, at least at high multipoles. They increase in size at smaller angular scales (i.e. 
large (.), since the APS becomes more sensitive to the 1-halo term. 


The left panel of Eig. 23 also demonstrates that the expected cross-correlation be¬ 
tween cross-correlation and the emission of unresolved blazars is more than 2 orders 
of magnitude smaller than the one with the other classes of astrophysical sources in 
Ref. [75], i.e. SEGs and MAGNs. As commented in the previous section, this is because 
the Fermi LAT has detected most of the blazars populating the volume probed by cosmic 
shear (i.e. z < 2 [75]) and, thus, the window functions Hbiazars(x) W^ix) do not 


L-y{M) relation enters in the computation of the 2-halo term only through the bias factor, 
bx{Y{M)), which is, generally, of 0{l). 
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Figure 23: Left: Cross-correlation APS between cosmic shear and gamma rays from blazars (red), 
MAGNs (pink), and SFGs (orange). Among the curves with the same color, different lines correspond to 
different choices for the relation: solid lines represent the fiducial models considered in Ref. [75], 

while dashed ones indicate extreme scenarios assumed to bracket the uncertainty on L^{M). Right: 
Gross-correlation APS between cosmic shear and DM-induced gamma-ray emission. Blue lines corre¬ 
spond to an annihilating DM candidate and the green one to a decaying DM candidate. The different 
blue lines represent different scenarios for low-mass DM (sub)halos (see Ref. [75] for details). The mass 
of the DM particle is taken to be 100 GeV (200 GeV) in the case of annihilating (decaying) DM. The 
annihilation cross section is fixed at 3 x 10“^®cm®s“^ and the decay lifetime at 3 x lO^^s. Annihilations 
and decays entirely into bb are assumed. Taken from Ref. [75]. 


overlap significantly. On the other hand, only a limited number of MAGNs and SFGs 
have been observed by the Fermi LAT (see also Secs. 2.2.2 and 2.2.3), so that the emis¬ 
sion of their unresolved counterparts is still characterized by a notable cross-correlation 
with the lensing signal. 

The right panel of Fig. 23 shows the expected cross-correlation APS with the gamma- 
ray emission produced by annihilating DM (blue lines) or decaying DM (green line). 
Different blue lines corresponds to different models for the description of low-mass DM 
halos. As seen in the previous sections, this uncertainty can have a significant impact on 
the properties of the DM-induced emission. Ref. [75] estimates that different description 
of the DM (sub)halos below the mass resolution of N-body simulations can lead to 
an uncertainty as large as a factor of 100 on the expected cross-correlation APS."^® In 
the most optimistic scenario considered in Ref. [75], the cross-correlation APS from 
annihilating DM is of the same order of that expected from astrophysical sources. This 
demonstrates how effective this observable can be for the detection of the DM component 
in the DGRB. 

Ref. [75] estimates that the measurement of the cross-correlation between the data 
of the Dark Energy Survey [439] and those of the Fermi LAT (after 5 years of operation) 
has the potential to detect an annihilating DM particle with an annihilation cross section 
smaller that the thermal value of 3 x 10“^®cm^s“^ for DM masses up to 300 GeV (for 


'^^Note, however, that the “NS” model (dotted blue line in Fig. 23) is probably underestimating the 
signal from DM since it assumes no contribution from DM subhalos and a quite large Mmin = 1O^M0. 
On the other hand, the “HIGH” scenario (dashed blue line) is based on power-law extrapolations from 
Ref. [276] and, thus, it likely overestimates the DM signal (see Sec. 2.3). 
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annihilations into b quarks). The result refers to a very optimistic subhalo boost and the 
predicted signal decreases by a factor of ~ 10 in the case of a more realistic description 
of DM subhalos (see right panel of Fig. 23. Prospects improve significantly with the 
inclusion of data from the forthcoming Euclid mission (expected for 2020) [440, 441]. 
Indeed, the cross correlation of Euclid data with those of a hypothetical future gamma- 
ray telescope with improved performances with respect to the Fermi LAT^^ has the 
potential to detect a DM component for DM masses up to the TeV scale (assuming a 
thermal annihilation cross section and an optimistic subhalo boost). 
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Figure 24: The cross-correlation signal of cosmic shear and the DGRB. The 4 panels corresponds to 
the 4 sky-patches W1-W4 covered by CFHTLenS. Red points show the results obtained using tangential 
shear (indicated as yt), while black points are for a component of the shear (yx in the legends) rotated by 
45° with respect to yt. The error bars indicate the standard deviation, estimated from 500 randomized 
shear catalogs. The quantifies the significance of the signal with respect to the statistical error. Taken 
from Ref. [71]. 


Such an optimal prospect for a DM detection is possible thanks to the “tomographic 


®°This improved version the Fermi LAT is described by a wide energy range, i.e. from 300 MeV to 1 
TeV. It is also assumed to achieve ~ 2.5 times more exposure (defined in Ref. [75] as the product of the 
effective area and the observation time) than 5 years of Fermi LAT data and to drastically improve the 
angular resolution to a value of 0.027° over the whole energy range. Its held of view is similar to that of 
the Fermi LAT. 
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spectral approach” employed in Ref. [75] , which combines spectral information with the 
study of the dependence of Cf’'^ on redshift (i.e. tomography). Such a technique provides 
an excellent sensitivity to DM-induced emission, even if the intensity or auto-correlation 
APS of such a component are only subdominant. 

The first (and, to date, only) measurement of the cross-correlation between DGRB 
and cosmic shear is performed in Ref. [71]. The authors make use of the data of 
CFHTLens from Ref. [442]. The survey detected more than 5 million galaxies in 4 
patches, covering an area of 154 square degrees. The corresponding shear signal is cor¬ 
related in Ref. [71] with 65 months of Fermi LAT data from the same region in the 
sky. Only photons between 1 and 500 GeV are considered in the analysis. A model 
for the diffuse Galactic foreground (determined in Ref. [71] from the gamma-ray data) 
is subtracted from the total gamma-ray emission and 2FGL sources are masked, before 
computing the cross-correlation. 

The estimator considered in Ref. [71] for the 2-point correlation function is 


{1 + 


(45) 


where nj is the number of gamma rays in the pixel centered on direction (pi, and ej 
is the shear-induced tangential ellipticity in pixel (pj = (pi + 'dj. The factors Wj and 
K{'&) depend on the precision in the estimation of ej (see Refs. [443, 71]) and the sum 
runs over all the pairs of pixels available. Fig. 24 shows the 2-point cross-correlation 
function (binned in 10 logarithmic bins with Alogj^Q?? = 0.2) for the 4 GFHTLenS sky 
patches. The measured 2-point correlation function using Eq. (45) are showed in red, 
while the black points are obtained from another shear component rotated by 45° with 
respect to the tangential one. In the case of a perfect measurement of the shape of the 
galaxies and of no intrinsic alignment, there should be no cross-correlation with this 
rotated data set. Thus, the black points in Fig. 24 represent a control sample with 
no cross-correlation. Note that black and red points in Fig. 24 are compatible with 
each other within errors, proving that no significant cross-correlation is present between 
GFHTLenS and Fermi LAT. This result can be translated into an upper limit on the 
DM annihilation cross section as a function of its mass. Ref. [71] excludes cross sections 
as low as the thermal value of 3 x 10“^®cm^s“^ for DM masses smaller than ~ 10 GeV 
(in the case of annihilation into and of the optimistic subhalo boost model of 

Ref. [276]). 


5.3. The cross-correlation with other tracers 

In addition to galaxy catalogs and cosmic shear, it is also possible to consider other 
observables that trace the LSS of the Universe. Below, we highlight some of these studies. 

Ref. [399] cross-correlates the DGRB inferred from the first 21 months of Fermi LAT 
data with the GMB measured by WMAP7 [444]. Such a measurement has the potential 
to probe the properties of dark energy through the detection of the so-called integrated 
Sachs-Wolfe effect [445]. This arises when the LSS gravitational potential changes with 
time during a cosmic era dominated by dark energy, as, e.g., in the local Universe. 
Additional anisotropies are induced in the GMB, which are expected to correlate with 
the LSS and, thus, potentially with the DGRB. In Ref. [399], the 2-point correlation 
function is computed, similarly to what was done for galaxy catalogs in Sec. 5.1. Their 
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results are compatible with a null detection, due to the large statistical errors.®^ Overall, 
Ref. [399] shows that the goal of detecting the integrated Sachs-Wolfe effect by cross- 
correlating the DGRB with the CMB is not unrealistic, but beyond the reach of the 
limited gamma-ray sample considered in Ref. [399]. 

More recently, Ref. [446] computed, for the first time, the cross-correlation of the 
DGRB with the so-called “leasing potential” of the GMB: the gravitational leasing 
induced by LSS imprints some distortions on the anisotropy pattern of the GMB, in 
such a way that the radiation detected by the Planck satellite [447] today is not exactly 
the one emitted at recombination. A statistical analysis of the non-Gaussianity of the 
GMB allows to reconstruct the leasing potential responsible for such perturbations [448, 
449, 450] . The first all-sky map of the CMB leasing potential has been recently reported 
by the Planck collaboration [451]. The signal is mainly contributed by structures at 
z ~ 2 and it exhibits an auto-correlation APS that peaks at £ ~ 20 — 30. 

Ref. [446] cross-correlates the sky-map of the CMB leasing potential with 68 months 
of Fermi LAT data, after having removed the diffuse Galactic emission. Six energy 
bins are considered, between 700 MeV and 300 GeV. The region at \b\ < 25° along the 
Galactic plane and a l°-circle around each source in the 2FGL are masked, together with 
the baseline 70% Galactic mask from Ref. [451]. The signal region is defined between 
^ = 40 and 400. A detection with a significance of 3.2cr is reported in the low-multipole 
region (£ < 160). No signal is present above i = 160 (see data points in Fig. 25). 
The signal at low multipoles is compared to the predicted cross-correlation between 
the lensing potential and the gamma-ray emission from 4 classes of unresolved sources, 
namely FSRQs, BL Lacs, SFGs and MAGNs. The LFs of these populations are fixed 
to the best-ht models from Refs. [22, 23, 29, 164, 33] (see Secs. 2.2.1, 2.2.2 and 2.2.3). 
The combined emission of these 4 source classes (solid black line in Fig. 25) reproduces 
fairly well the experimental data. More precisely, BL Lacs, SFGs and MAGNs (red, 
orange and green lines in Fig. 25, respectively) contribute more or less in equal parts to 
the cross-correlation signal, while FSRQs (blue line) are subdominant. This model also 
provides a good fit to the DGRB energy spectrum and auto-correlation APS. 


6. The science of the Diffuse Gamma-Ray Background 

In this review we have summarized the current knowledge on the Diffuse Gamma- 
Ray Background (DGRB). The DGRB is what remains in the gamma-ray sky after 
the subtraction of the diffuse Galactic foreground and of the resolved sources. It is 
interpreted as the cumulative emission of the objects that are not bright enough to be 
detected individually. 

Since its first detection in 1972 by the OSO-3 satellite [1], this emission has been 
deeply investigated in an attempt to understand its composition. The Fermi LAT satel¬ 
lite, in operation since 2008, has greatly improved our understanding of the DGRB. This 
emission is measured as an isotropic template in a multi-component fit to the Fermi LAT 
data. The fit also includes a model for the resolved sources and for the diffuse Galactic 
foreground. The most recent measurement of the DGRB energy spectrum is reported 
in Ref. [9] and it covers almost 4 orders of magnitude in energy, from 100 MeV to 820 
GeV. Our imperfect knowledge of the Galactic foregrounds represents the main source 


®^The lack of a significant detection is compatible with the expectation if the DGRB is composed by 
unresolved sonrces (parametrized as in Ref. [399]). 
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Figure 25: Cross-correlation APS CJ'^ between the DGRB and the CMB lensing potential, as a 
function of the multipole for gamma-ray energies E > 1 GeV. The measurement is averaged (linearly 
in terms of in multipole bins of t\l = 60, starting at £ = 40. Points report the minimum- 

variance combination of the measurement in the individual energy bins (assuming a spectrum oc 
as described in Ref. [446]. Four different analyses are shown. They arise from the combination of two 
lensing maps (from the 2013 and 2015 release of Planck data) and two gamma-ray point-source masks 
(masking sources in 2FGL or 3FGL). The benchmark theoretical model, shown in black, is the sum of the 
contributions from BL Lacs (red), FSRQs (blue), MAGNs (magenta) and SFGs (orange), multiplied by a 
normalization factor A”'’' = 1.35. Two generic models (labeled “GO.l” and “G2”) with Gaussian window 
functions are also shown. The peak of the Gaussian is at zo = 0.1 {zo = 2), with a dispersion = 0.1 
{az = 0.5) for GO.l (G2). In the inset, the intensity energy spectrum is shown for the Fermi LAT 
measurement (black data point, labeled “EGB”) and for the model predictions. Taken from Ref. [446]. 


of uncertainty in the analysis and it induces a systematic error on the DGRB intensity 
of ~ 15 — 30%, depending on the energy range considered. 

Before the Fermi LAT, the energy spectrum was the only source of information 
available on the DGRB. However, the scenario drastically changed in 2012 when, for 
the first time, the Fermi LAT measured also the angular power spectrum (APS) of 
anisotropies in the DGRB [66]. The emission was found to exhibit a Poissonian APS 
in the multipole range between i = 155 and 504, with a significance as large as 7.2a 
between 1.99 and 5.0 GeV. This signal is independent of energy between 1 and 50 GeV. 

Recently, cross-correlations of the DGRB with different data sets have also revealed to 
be a powerful tool to unveil the composition of the DGRB. Ref. [70] reported a significant 
cross-correlation of the DGRB with 4 out of the 5 galaxy catalogs considered, namely 
the optically-selected quasares of the Sloan Digital Sky Survey (SDSS) Data Release 
6 in Ref. [425], the IR galaxies of the 2 Micron All-Sky Survey (2MASS) Extended 
Source Catalog from Ref. [426], the main galaxy sample of SDSS Data Release 8 from 
Ref. [429] and the radio sources from the NRAO VLA Sky Survey (NVSS) [427]. These 
cross-correlation signals, obtained after the analysis of 60 months of Fermi LAT data, 
are localized at small angles (below few degrees) with significances that range between 
3 (t (in the case of SDSS Data Release 8 main galaxies) and more than lOcr (for the 


68 












































NVSS catalog), for gamma-ray energies above 1 GeV. The cross-correlation with NVSS, 
however, is most likely contaminated by a 1-halo term not related to the Large Scale 
Structure (LSS) of the Universe. 

In Ref. [71], the authors measured, for the first time, the cross-correlation of the 
DGRB with the cosmic shear induced by the gravitational lensing detected in the 
Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS) [442]. Their results are 
compatible with a null cross-correlation signal in the angular range between 1 and 100 
arcmin. Ref. [71] proves, though, that measuring the cross-correlation with the cosmic 
shear signal is not an unrealistic goal. Indeed, more promising results are expected from 
the cross-correlation between the DGRB and the data expected from the Dark Energy 
Survey (DES) or from Euclid [75]. Finally, Ref. [446] reported a 3.2cr detection of the 
cross-correlation between the DGRB and the lensing potential of the CMB measured by 
the Planck Gollaboration [451]. The signal is localized in the multipole region between 
£ = 40 and 160, for gamma-ray energies between 700 MeV and 300 GeV. 

The increased amount of observational data on the DGRB has allowed significant 
progress in the modeling of its composition. The DGRB is interpreted as the cumu¬ 
lative emission of unresolved sources, e.g. blazars, misaligned Active Galactic Nuclei 
(MAGNs), star-forming galaxies (SFGs) and MilliSecond Pulsars (MSPs). It, therefore, 
represents a reservoir of invaluable information on these astrophysical sources as it may 
be the only way to study the emission of objects that are too faint to be detected individ¬ 
ually. In particular, the DGRB can potentially determine the faint end of the luminosity 
function of the aforementioned populations, i.e. a goal that would be quite difficult, if 
not impossible, to achieve otherwise. 

Before the Fermi LAT, the predictions for the contribution of unresolved blazars to 
the DGRB were affected by large uncertainties [94, 10, 11, 12, 13, 15, 30, 16]. Nowadays, 
the wealth of new information on the DGRB, combined with the population studies of 
resolved blazars performed by the Fermi LAT Gollaboration, have established that this 
source class cannot account for more than (20 ± 4)% of the DGRB in the energy range 
between 0.1 and 100 GeV [25]. The subclass of blazars responsible for the bulk of the 
blazar contribution varies depending on the energies considered. Indeed, Refs. [23, 24, 25] 
showed that unresolved high-synchrotron-peaked BL Lacs can explain the whole DGRB 
at energies above ~ 100 GeV. At lower energies, astrophysical populations other than 
blazars are required. Yet, the modeling of these other gamma-ray emitters, such as 
namely MAGNs, SFGs and MSPs, is not as robust as that of blazars. This is caused by 
difficulty in performing reliable population studies with the limited sample of resolved 
sources currently available in the gamma-ray range. Generally, it is useful to assume 
a correlation between luminosities at different wavelengths (i.e. gamma-ray with radio 
frequencies [28, 29] in the case of MAGNs, or gamma rays and infra-red light [33] for 
SFGs). 

The overall picture indicates that the 4 classes of astrophysical sources mentioned 
above are enough to explain the totality of the DGRB energy spectrum measured in 
Ref. [9] (see Fig. 9). As for the APS of DGRB anisotropies. Refs. [67, 68, 63] prove that 
unresolved blazars alone (more specifically, high-synchrotron-peak BL Lacs) can account 
for the whole APS reported in Ref. [66]. These two important results can be reconciled 
by the fact that the blazar component is produced by a relatively small number of 
bright but unresolved objects and, thus, it gives rise to significant anisotropies. On the 
contrary, the other source classes produce fairly isotropic cumulative emission, as their 
members are more numerous and fainter. 

The picture gets more complicated, though, when considering also the measured 
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cross-correlation with LSS tracers: the sum of unresolved blazars, MAGNs and SFGs 
provides a good fit to the cross-correlation APS detected between the DGRB and the 
CMB lensing potential [446]. It is also compatible with the lack of significant cross¬ 
correlation with the GFHTLenS cosmic shear [71]. However, Ref. [430] finds that the 
model that best fits the two-point correlation functions measured in Ref. [70] with 5 
galaxy catalogs can only account for ~ 30% of the DGRB intensity reported in Ref. [9]. 
It is still unclear if such a limitation of the astrophysical interpretation of the DGRB 
can be alleviated by a more sophisticated modeling of its components. Also, it will be 
interesting to verify whether a similar scenario will be confirmed by new surveys or by 
more complete data releases of the catalogs currently employed. Alternatively, one will 
be forced to supplement the model with another component, which does not correlate 
with the LSS, such as gamma-ray emission associated with the DM halo of the MW or 
with its DM substructures. 

Quantifying the contribution of known astrophysical populations automatically con¬ 
strains the intensity of other potential contributors to the DGRB. We briefly presented, 
among others, the case of clusters of galaxies. Type la supernovae and of Ultra-High- 
Energy Cosmic Rays interacting with background radiation. However, the most studied 
scenario is that of a potential gamma-ray emission from Dark Matter (DM) annihilation 
or decay. Since no DM signal has been undoubtedly detected in the gamma-ray sky so 
far, it is expected that most of this hypothetical gamma-ray emission will contribute to 
the DGRB. Searching for DM in the DGRB has the advantage that the DM-induced 
component is sourced by the emission coming from all DM halos and subhalos around 
us. It will, thus, depend on the ensemble-averaged properties of the DM halo population, 
which can be inferred from Wbody cosmological simulations [452, 343, 267, 453, 454] 
or predicted by the theory of structure formation [200, 201]. This is intrinsically dif¬ 
ferent from the observation of a specific target, e.g. the Galactic Center or a dwarf 
Spheroidal satellite galaxy. Indeed, each of these targets could be very peculiar and 
deviate considerably from the ensemble average, potentially hamper the interpretation 
of any data. Another benefit of using the DGRB to search for DM is that a potential 
DM signal contributing to the DGRB would be sensitive to the process of assembly of 
DM halos and their subsequent evolution. This kind of information would be difficult 
(if not impossible) to extract by observing individual targets in the sky. In this regard, 
the DGRB may be the only cosmological non-gravitational probe of DM. In addition, 
it represents a fundamental source of complementary information in the study of any 
claimed DM signal. 

The intensity of this DM-induced cosmological emission depends on the properties 
of the DM particle, e.g. its mass, annihilation cross section or decay lifetime. Also, it 
rests on the abundance and properties of DM structures. Our understanding of DM 
(sub)halos heavily relies on the results of Al-body cosmological simulations. Yet, as 
of today, even the simulations with the highest resolution are far from resolving the 
whole DM halo hierarchy down to the predicted Mmin. The properties of low-mass 
DM structures, thus, need to be inferred by extrapolating the characteristics of their 
more massive counterparts, that are well resolved in the simulations. Heuristic power- 
law extrapolations, which predict very bright low-mass DM halos, have been commonly 
used in the literature [274, 56, 275, 276, 339, 217]. However, recent high-resolution 
simulations of the smallest DM halos [277, 279] suggest that those extrapolations are 
not well motivated, favoring, instead, models that yield a more moderate DM-induced 
gamma-ray emission from low-mass structures [267, 280, 281]. Overall, this results in 
a DM-induced cosmological signal which is substantially weaker than the one obtained 
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when assuming power-law extrapolations. Also, and perhaps even more importantly, 
the improved knowledge provided by Refs. [267, 280, 281] on the structural properties of 
the smallest DM halos has considerably reduced the theoretical uncertainty associated 
with the DM contribution to the DGRB down to a factor of ~ 20 [64]. 

Not unexpectedly, the measured DGRB energy spectrum can be used to constrain 
the intensity of the DM-induced emission and, thus, to derive upper (lower) limits on 
the annihilation cross section (decay lifetime). Figs. 11 and 12 summarize some of 
these results. For annihilating DM, the so-called “sensitivity-reach” limits derived in 
Ref. [64] from their hducial Halo Model exclude thermal annihilation cross sections of 
3 X 10“^®cm^s“^ for masses below ~ 100 GeV in the case of annihilations into b quarks. 
When compared to other indirect searches for DM, the upper limits inferred from the 
DGRB represent the strongest constraints on (crv) currently available, for DM masses 
below ~ 1 TeV. A more conservative statistical analysis would increase the upper limits 
by a factor of ~ 10.^^ In the case of decaying DM, the lower limits on r obtained in 
Ref. [386] from the DGRB measurement in Ref. [8] exclude decay lifetimes smaller than 
~ 3 X 10^’^ s for < 2 TeV and decays into At larger masses, the strongest 

lower limit comes from Ref. [387], which excludes lifetimes as large as 10^® s. When 
compared to other searches of decaying DM, these limits represent the most constraining 
information available on r, up to DM masses of 20 TeV. 

A similar strategy can be used when comparing the APS measured in Ref. [66] to the 
predicted anisotropies in the DM-induced emission. Nevertheless, the latter turns out 
to be quite isotropic [56, 62, 339] and, thus, the corresponding upper limits, although 
competitive with other indirect DM searches, are less stringent than those inferred from 
the DGRB [339, 420]. The cross-correlation with galaxy catalogs and with cosmic shear 
are also very promising strategies to constrain (or even to detect) a potential DM signal 
in the DGRB [73, 74, 72, 75]. This is possible mainly thanks to the fact that both cross¬ 
correlations are particularly sensitive to the way the matter is distributed in the local 
Universe and, in particular, to the most massive DM halos. We remind that unresolved 
blazars, the main contributors to the auto-correlation APS, do not populate neither the 
local volume nor the largest DM halo masses. A model of the DGRB including both 
astrophysical sources and annihilating DM is used in Ref. [430] to describe measurement 
of the 2-point correlation function reported in Ref. [70]. The analysis excludes DM 
candidates with annihilation cross sections larger than the thermal value for masses 
below 40 GeV, for annihilations into b quarks and a moderate value of the subhalo boost. 
When assuming the same value of the subhalo boost, the obtained upper limit is currently 
the strongest one among all those derived from DGRB data, including measurements of 
the DGRB intensity and of the auto-correlation APS. 

The energy spectrum, auto-correlation and cross-correlation APS of the DGRB are 
sensitive to different characteristics of the sources contributing to the emission. Consid¬ 
ering at the same time all three observables provide a very powerful handle to reconstruct 
the composition of the DGRB. Needless to say, such an ambitious goal requires ample 
data sets, including information from wavelengths other than the gamma-ray energy 
range. Fortunately, a great wealth of new observational information is expected in the 
near future: the Fermi LAT will continue gathering data until, at least, 2016. The 
Cherenkov Telescope Array (CTA), expected to be in operation by 2020, will improve 
by a factor ~ 10 the sensitivity of current Cherenkov telescopes in the energy range be- 


®^Note also that both the fiducial sensitivity-reach limits and the more conservative ones in Ref. [64] 
are affected by an uncertainty of a factor of ~ 3. 
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tween a few dozens of GeV and a few dozens of TeV. Its complementarity with respect to 
the Fermi LAT will allow an improved precision in the determination of the DGRB en¬ 
ergy spectrum in the sub-TeV energy range, as well as the extension of the measurement 
beyond the TeV. In addition, GTA will perform the first survey of a significant portion 
of the sky at these very-high energies [455, 456]. Combined with the Fermi LAT data 
gathered since Ref. [66], it will be possible to extend the data on the auto-correlation 
APS to higher energies and to decrease the size of the bins in energy. On the other 
hand, the cross-correlation of the DGRB with LSS tracers will hugely benefit for the 
imminent release of the data gathered by the DES [439] during its first year of operation. 
The near future will also see the advent of the next generation of galaxy catalogs, e.g. 
the extended Baryon Oscillation Spectroscopy Survey (eBOSS)^^ and the Dark Energy 
Spectroscopic Instrument (DESI, formerly BigBoSS) [457]. On a longer scale, Euclid 
will offer weak lensing measurements with unprecedented precision after 2020 [440]. 

The increased observational data available on the DGRB will also be accompanied 
by significant progress in the modeling of its contributors. Improving our understanding 
of the emission of blazars, SEGs and MAGNs will alleviate the degeneracies currently 
affecting our interpretation and will reduce the uncertainty associated with each com¬ 
ponent. A fully multi-wavelength approach is required to achieve such a goal: in the 
X-ray band the Nuclear Spectroscopic Telescope Array (NuSTAR) [458] has been re¬ 
cently launched and ASTRO-H [459, 460] will follow within this year (2015). Infra-red 
data from Herschel and the Wide-field Infrared Survey Explorer (WISE)^^ are already 
public and the James Webb Space Telescope (JWST)^^ is expected to be launched in 
2018. Regular observation in radio has started with Low-Erequency Array (LOEAR) 
a pathfinder for the Square Kilometer Array (SKA) planned for 2020. 

We end the review by noting that the IceGube Gollaboration has recently reported the 
detection of the first extraterrestrial neutrinos [461, 462, 463]. The 37 observed neutrino 
events represent an excess over the atmospheric background extending to the PeV scale, 
with a significance of more than ha [463]. Even if the origin of these neutrinos is not 
clearly established yet (see Refs. [464, 465] and references therein), a possibility is that 
they originate from a diffuse neutrino flux similar to the DGRB.^® If this interpretation 
was confirmed, many of the techniques discussed in this review could also be adopted 
to investigate this diffuse neutrino flux. Indeed, many of the sources contributing to 
the DGRB are expected to emit neutrinos as well [466, 467, 468, 469]. This implies 
that any constraint on their neutrino emission would indirectly constrain also their 
contribution to the DGRB (and viceversa). The development of a fully multi-messenger 
approach is certainly a very tantalizing possibility that has been already considered, e.g., 
in Refs. [470, 471, 472, 207, 473, 161]. 

In conclusion, the DGRB is a fundamental component of the gamma-ray sky, whose 
exact composition still remains unveiled. The recent rapid growth of available data 
on the DGRB (mostly thanks to the outstanding performance of the Fermi LAT) has 
triggered an increased attention from the scientihc community. Astrophysicists and 
astroparticle physicists aim at reconstructing the composition of the DGRB to infer 


®®https: //www.sdss3.org/future/eboss.php 
®^http: / / wise.ssl.berkeley.edu/ 
http://www .jwst.nasa.gov/ 

®®http: / / www.lofar.org/ 

®^http: / / www.skatelescope.org/ 

®®To stress the similarity with the DGRB, we propose to call this neutrino emission, the Diffuse 
Neutrino Background (DNB). 
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novel information on the sources contributing to the emission, especially in their low- 
luminosity regime. New data sets are already available (or will be soon) which can 
provide a significant progress to the common goal of dissecting the true nature of the 
DGRB. By summarizing where we stand on our current understanding of this emission, 
with this review we hope to have offered a useful reference to those who will analyze and 
interpret the data to come, as well as to help finding new avenues and opportunities for 
further research. 
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